From 19fb9ea598a236bffad758c879d9040ea22aabc2 Mon Sep 17 00:00:00 2001 From: g_giud Date: Sat, 27 May 2017 21:15:31 -0400 Subject: [PATCH 01/54] partial currents in python interface Added cell_from and cell_to tallies --- openmc/filter.py | 73 +++++++++++++++++++++++++++++++++++++++++++++-- openmc/tallies.py | 21 +++++++++++++- 2 files changed, 91 insertions(+), 3 deletions(-) diff --git a/openmc/filter.py b/openmc/filter.py index 676bc66e17..063961d814 100644 --- a/openmc/filter.py +++ b/openmc/filter.py @@ -15,7 +15,7 @@ import openmc.checkvalue as cv _FILTER_TYPES = ['universe', 'material', 'cell', 'cellborn', 'surface', 'mesh', 'energy', 'energyout', 'mu', 'polar', 'azimuthal', - 'distribcell', 'delayedgroup', 'energyfunction'] + 'distribcell', 'delayedgroup', 'energyfunction', 'celltocell', 'cellfrom', 'cellto'] _CURRENT_NAMES = {1: 'x-min out', 2: 'x-min in', 3: 'x-max out', 4: 'x-max in', @@ -335,6 +335,9 @@ class Filter(object): """ if filter_bin not in self.bins: + print('Failed to find') + print(filter_bin) + print(self.bins) msg = 'Unable to get the bin index for Filter since "{0}" ' \ 'is not one of the bins'.format(filter_bin) raise ValueError(msg) @@ -533,7 +536,66 @@ class CellFilter(IntegralFilter): """ +class CellToCellFilter(IntegralFilter): + """Bins tally on which couple of cells the neutrons went from and to. /CHANGE/ + Parameters + ---------- + bins : Iterable of Integral + openmc.Cell IDs. + + Attributes + ---------- + bins : Iterable of Integral + openmc.Cell IDs. + num_bins : Integral + The number of filter bins + stride : Integral + The number of filter, nuclide and score bins within each of this + filter's bins. + + """ + +class CellFromFilter(IntegralFilter): + """Bins tally on which couple of cells the neutrons came from. /CHANGE/ + + Parameters + ---------- + bins : Integral or Iterable of Integral + openmc.Cell IDs. + + Attributes + ---------- + bins : Integral or Iterable of Integral + openmc.Cell IDs. + num_bins : Integral + The number of filter bins + stride : Integral + The number of filter, nuclide and score bins within each of this + filter's bins. + + """ + +class CellToFilter(IntegralFilter): + """Bins tally on which couple of cells the neutrons went to. /CHANGE/ + + Parameters + ---------- + bins : Integral or Iterable of Integral + openmc.Cell IDs. + + Attributes + ---------- + bins : Integral or Iterable of Integral + openmc.Cell IDs. + num_bins : Integral + The number of filter bins + stride : Integral + The number of filter, nuclide and score bins within each of this + filter's bins. + + """ + class CellbornFilter(IntegralFilter): """Bins tally events based on the cell that the particle was born in. @@ -901,8 +963,12 @@ class RealFilter(Filter): return np.allclose(self.bins, other.bins) def get_bin_index(self, filter_bin): - i = np.where(self.bins == filter_bin[1])[0] + i = np.where(abs(self.bins -filter_bin[1]) < 1e-10)[0] if len(i) == 0: + print('Tally bins', self.bins) + print('Bins searched for:', filter_bin) + print('Search result', i) + print(self.bins - filter_bin[1]) msg = 'Unable to get the bin index for Filter since "{0}" ' \ 'is not one of the bins'.format(filter_bin) raise ValueError(msg) @@ -940,6 +1006,9 @@ class EnergyFilter(RealFilter): def get_bin_index(self, filter_bin): # Use lower energy bound to find index for RealFilters + + #print(self.bins, filter_bin) + deltas = np.abs(self.bins - filter_bin[1]) / filter_bin[1] min_delta = np.min(deltas) if min_delta < 1E-3: diff --git a/openmc/tallies.py b/openmc/tallies.py index 0bf87de101..f3af788a65 100644 --- a/openmc/tallies.py +++ b/openmc/tallies.py @@ -239,6 +239,9 @@ class Tally(object): num_bins = 1 for self_filter in self.filters: + + #print("Looping through filters, num_bin", self_filter.num_bins, num_bins) + num_bins *= self_filter.num_bins return num_bins @@ -1215,6 +1218,7 @@ class Tally(object): # Find the equivalent Filter in this Tally's list of Filters filter_found = self.find_filter(filter_type) + #print('searched for filter', filter_type, 'and found', filter_found) # Get the index for the requested bin from the Filter and return it filter_index = filter_found.get_bin_index(filter_bin) @@ -1376,6 +1380,12 @@ class Tally(object): else: bins = self_filter.bins + #print('in get_filter_indices') + #print('arguments', filters, filter_bins) + #print('params') + #print(bins) + #print(type(self_filter), self_filter) + # Initialize a NumPy array for the Filter bin indices filter_indices.append(np.zeros(len(bins), dtype=np.int)) @@ -1520,6 +1530,8 @@ class Tally(object): e.g., if the score(s) do not match those in the Tally. """ + #print("\n In get_values") + #print("Tally shape", self.shape) # Ensure that the tally has data if (value == 'mean' and self.mean is None) or \ @@ -1529,6 +1541,13 @@ class Tally(object): (value == 'sum_sq' and self.sum_sq is None): msg = 'The Tally ID="{0}" has no data to return'.format(self.id) raise ValueError(msg) + + #print('\n arguments') + #print(filters) + #print(filter_bins) + #print('\n Tally characteristics') + #print(self) + #print(self.name, self.filters, self.scores,self.num_filter_bins) # Get filter, nuclide and score indices filter_indices = self.get_filter_indices(filters, filter_bins) @@ -1537,7 +1556,7 @@ class Tally(object): # Construct outer product of all three index types with each other indices = np.ix_(filter_indices, nuclide_indices, score_indices) - + # Return the desired result from Tally if value == 'mean': data = self.mean[indices] From ff6f05c4b7431e484bedf63b897fec3ecfeca46d Mon Sep 17 00:00:00 2001 From: g_giud Date: Sat, 27 May 2017 21:18:03 -0400 Subject: [PATCH 02/54] partial current tallies in fortran src coded when they should be triggered, when doing surface crossings, and in which form , as cell_from and cell_to filters --- src/constants.F90 | 13 +- src/geometry.F90 | 19 ++- src/global.F90 | 4 +- src/input_xml.F90 | 67 ++++++++- src/particle_header.F90 | 13 +- src/tally.F90 | 141 +++++++++++++++++- src/tally_filter.F90 | 315 +++++++++++++++++++++++++++++++++++++++ src/tally_initialize.F90 | 1 + src/tracking.F90 | 11 ++ 9 files changed, 573 insertions(+), 11 deletions(-) diff --git a/src/constants.F90 b/src/constants.F90 index c44896b316..b5ae55796a 100644 --- a/src/constants.F90 +++ b/src/constants.F90 @@ -288,7 +288,8 @@ module constants ! Tally type integer, parameter :: & TALLY_VOLUME = 1, & - TALLY_SURFACE_CURRENT = 2 + TALLY_SURFACE_CURRENT = 2, & + TALLY_CELL_TO_CELL = 3 ! Tally estimator types integer, parameter :: & @@ -304,7 +305,7 @@ module constants EVENT_ABSORB = 2 ! Tally score type - integer, parameter :: N_SCORE_TYPES = 24 + integer, parameter :: N_SCORE_TYPES = 25 integer, parameter :: & SCORE_FLUX = -1, & ! flux SCORE_TOTAL = -2, & ! total reaction rate @@ -329,7 +330,8 @@ 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_CELL_TO_CELL_TYPE = -25 ! cell to cell partial current ! Maximum scattering order supported integer, parameter :: MAX_ANG_ORDER = 10 @@ -352,7 +354,7 @@ module constants integer, parameter :: NO_BIN_FOUND = -1 ! Tally filter and map types - integer, parameter :: N_FILTER_TYPES = 14 + integer, parameter :: N_FILTER_TYPES = 15 integer, parameter :: & FILTER_UNIVERSE = 1, & FILTER_MATERIAL = 2, & @@ -367,7 +369,8 @@ module constants FILTER_POLAR = 11, & FILTER_AZIMUTHAL = 12, & FILTER_DELAYEDGROUP = 13, & - FILTER_ENERGYFUNCTION = 14 + FILTER_ENERGYFUNCTION = 14, & + FILTER_CELL_TO_CELL = 15 ! Mesh types integer, parameter :: & diff --git a/src/geometry.F90 b/src/geometry.F90 index 4856159913..cb9236ba2b 100644 --- a/src/geometry.F90 +++ b/src/geometry.F90 @@ -11,7 +11,8 @@ module geometry use surface_header use stl_vector, only: VectorInt use string, only: to_str - use tally, only: score_surface_current + use tally, only: score_surface_current, & + &score_cell_to_cell implicit none @@ -396,6 +397,8 @@ contains integer :: i_surface ! index in surfaces logical :: found ! particle found in universe? class(Surface), pointer :: surf + + !print *, "cross surface called" i_surface = abs(p % surface) surf => surfaces(i_surface)%obj @@ -565,6 +568,11 @@ contains ! cells on the positive side call find_cell(p, found, surf%neighbor_pos) + + ! /CHANGE/ Score cell_to_cell +! print *, "Positive side, calling score_cell_to_cell" + call score_cell_to_cell(p, last_cell) + if (found) return elseif (p % surface < 0 .and. allocated(surf%neighbor_neg)) then @@ -572,6 +580,11 @@ contains ! cells on the negative side call find_cell(p, found, surf%neighbor_neg) + + ! /CHANGE/ Score cell_to_cell +! print *, "Negative side, calling score_cell_to_cell" + call score_cell_to_cell(p, last_cell) + if (found) return end if @@ -604,6 +617,10 @@ contains return end if end if + + ! /CHANGE/ Score cell_to_cell +! print *, "Searched all cells then calling score_cell_to_cell" + call score_cell_to_cell(p, last_cell) end subroutine cross_surface diff --git a/src/global.F90 b/src/global.F90 index 780f323f37..7b13dc4be5 100644 --- a/src/global.F90 +++ b/src/global.F90 @@ -155,10 +155,11 @@ module global type(SetInt) :: active_tracklength_tallies type(SetInt) :: active_current_tallies type(SetInt) :: active_collision_tallies + type(SetInt) :: active_cell_to_cell_tallies type(SetInt) :: active_tallies !$omp threadprivate(active_analog_tallies, active_tracklength_tallies, & !$omp& active_current_tallies, active_collision_tallies, & -!$omp& active_tallies) +!$omp& active_tallies, active_cell_to_cell_tallies) ! Global tallies ! 1) collision estimate of k-eff @@ -515,6 +516,7 @@ contains call active_tracklength_tallies % clear() call active_current_tallies % clear() call active_collision_tallies % clear() + call active_cell_to_cell_tallies % clear() call active_tallies % clear() ! Deallocate track_identifiers diff --git a/src/input_xml.F90 b/src/input_xml.F90 index a5c894d10c..03e7abc5ae 100644 --- a/src/input_xml.F90 +++ b/src/input_xml.F90 @@ -2989,6 +2989,8 @@ contains ! Get pointer list to XML and get number of filters call get_node_list(node_tal, "filter", node_filt_list) n_filters = size(node_filt_list) + +! print *, "Number of filters", n_filters ! Allocate filters array allocate(t % filters(n_filters)) @@ -3012,8 +3014,9 @@ contains end if n_words = node_word_count(node_filt, "bins") case ("mesh", "universe", "material", "cell", "distribcell", & - "cellborn", "surface", "delayedgroup") - if (.not. check_for_node(node_filt, "bins")) then + "cellborn", "surface", "delayedgroup", "celltocell",& + "cellfrom","cellto") ! /CHANGE/ + if (.not. check_for_node(node_filt, "bins")) then call fatal_error("Bins not set in filter on tally " & // trim(to_str(t % id))) end if @@ -3048,6 +3051,63 @@ contains end select ! Set the filter index in the tally find_filter array t % find_filter(FILTER_CELL) = j + + case ('celltocell') ! /CHANGE/ + ! Allocate and declare the filter type + allocate(CellToCellFilter::t % filters(j) % obj) + select type (filt => t % filters(j) % obj) + type is (CellToCellFilter) + ! Allocate and store bins + filt % n_bins = n_words / 2 ! 2 cells in each bin +! print *, "Filter n bins", n_words / 2 + allocate(filt % cells(n_words)) ! cell storage OR tally storage ??? + call get_node_array(node_filt, "bins", filt % cells) + ! node_filt is a pointer to filter + ! + end select + ! Set the filter index in the tally find_filter array + t % find_filter(FILTER_CELL_TO_CELL) = j + + ! Set estimator + t % estimator = ESTIMATOR_ANALOG + + case ('cellfrom') ! /CHANGE/ + ! Allocate and declare the filter type + allocate(CellFromFilter::t % filters(j) % obj) + select type (filt => t % filters(j) % obj) + type is (CellFromFilter) + ! Allocate and store bins + filt % n_bins = n_words +! print *, "Filter n bins", n_words + allocate(filt % cells(n_words)) + call get_node_array(node_filt, "bins", filt % cells) + ! node_filt is a pointer to filter + ! + end select + ! Set the filter index in the tally find_filter array + t % find_filter(FILTER_CELL_TO_CELL) = j + + ! Set estimator + t % estimator = ESTIMATOR_ANALOG + + case ('cellto') ! /CHANGE/ + ! Allocate and declare the filter type + allocate(CellToFilter::t % filters(j) % obj) + select type (filt => t % filters(j) % obj) + type is (CellToFilter) + ! Allocate and store bins + filt % n_bins = n_words +! print *, "Filter n bins", n_words + allocate(filt % cells(n_words)) + call get_node_array(node_filt, "bins", filt % cells) + ! node_filt is a pointer to filter + ! + end select + ! Set the filter index in the tally find_filter array + t % find_filter(FILTER_CELL_TO_CELL) = j + + ! Set estimator + t % estimator = ESTIMATOR_ANALOG case ('cellborn') ! Allocate and declare the filter type @@ -3783,6 +3843,9 @@ contains case ('events') t % score_bins(j) = SCORE_EVENTS + case ('partial_current') ! /CHANGE/ + t % type = TALLY_CELL_TO_CELL + t % score_bins(j) = SCORE_CELL_TO_CELL_TYPE case ('elastic', '(n,elastic)') t % score_bins(j) = ELASTIC diff --git a/src/particle_header.F90 b/src/particle_header.F90 index a7309393b2..c6d478540d 100644 --- a/src/particle_header.F90 +++ b/src/particle_header.F90 @@ -15,6 +15,9 @@ module particle_header !=============================================================================== type LocalCoord + + integer :: last_cell = NONE + ! Indices in various arrays for this level integer :: cell = NONE integer :: universe = NONE @@ -57,7 +60,13 @@ module particle_header real(8) :: wgt ! particle weight real(8) :: mu ! angle of scatter logical :: alive ! is particle alive? - + + ! Previous cell before crossing surface + integer :: last_cell + ! Particle coordinates + integer :: last_n_coord ! number of current coordinates +! type(LocalCoord) :: last_coord(MAX_COORD) ! coordinates for all levels + ! Pre-collision physical data real(8) :: last_xyz_current(3) ! coordinates of the last collision or ! reflective/periodic surface crossing @@ -172,6 +181,8 @@ contains elemental subroutine reset_coord(this) class(LocalCoord), intent(inout) :: this + this % last_cell = NONE + this % cell = NONE this % universe = NONE this % lattice = NONE diff --git a/src/tally.F90 b/src/tally.F90 index 248fb16642..dffd798f20 100644 --- a/src/tally.F90 +++ b/src/tally.F90 @@ -116,9 +116,12 @@ contains !######################################################################### ! Determine appropirate scoring value. + +! print *, "Score bin", score_bin select case(score_bin) - + + ! use default case for cell_to_from tallies, with t % estimator == ESTIMATOR_ANALOG case (SCORE_FLUX, SCORE_FLUX_YN) if (t % estimator == ESTIMATOR_ANALOG) then @@ -1177,7 +1180,16 @@ contains end if end if + case(SCORE_CELL_TO_CELL_TYPE) + ! Not using default because score_bin and event_MT thing + ! score_bin is -25, and event_MT seems to match the cell id + ! print*, p % event_MT, score_bin +! print *, "Scoring", p % last_wgt, flux + score = p % last_wgt * flux + case default + !print *, "Found right score type" + !print * , t % estimator if (t % estimator == ESTIMATOR_ANALOG) then ! Any other score is assumed to be a MT number. Thus, we just need ! to check if it matches the MT number of the event @@ -2079,6 +2091,13 @@ contains case (SCORE_EVENTS) ! Simply count number of scoring events score = ONE + + case(SCORE_CELL_TO_CELL_TYPE) + ! Not using default because score_bin and event_MT thing + ! score_bin is -25, and event_MT seems to match the cell id + ! print*, p % event_MT, score_bin +! print *, "Scoring", p % last_wgt, flux + score = p % last_wgt * flux end select @@ -2111,6 +2130,8 @@ contains integer :: num_nm ! Number of N,M orders in harmonic integer :: n ! Moment loop index real(8) :: uvw(3) + +! print *, "Expand and score, score bin", score_bin, "score", score, "filter_index", filter_index select case(score_bin) case (SCORE_SCATTER_N, SCORE_NU_SCATTER_N) @@ -2194,6 +2215,7 @@ contains case default +! print *, "Adding to total in expand-and-score", score !$omp atomic t % results(RESULT_VALUE, score_index, filter_index) = & t % results(RESULT_VALUE, score_index, filter_index) + score @@ -3019,6 +3041,120 @@ contains end subroutine score_collision_tally +!=============================================================================== +! SCORE_CELL_TO_CELL tallies partial currents from cell_to to cell_from +!=============================================================================== + + subroutine score_cell_to_cell(p, last_cell) + + type(Particle), intent(in) :: p + integer, intent(in) :: last_cell + + ! if particle has leaked + ! DF could preserve leakage aiight + !if p % alive = .false. then + +!end if + integer :: i + integer :: i_tally + integer :: i_filt + integer :: k ! loop index for nuclide bins + ! position during the loop + integer :: filter_index ! single index for single bin + integer :: i_nuclide ! index in nuclides array + real(8) :: filter_weight ! combined weight of all filters + real(8) :: dummy_real8 + type(TallyObject), pointer :: t + + ! A loop over all tallies is necessary because we need to simultaneously + ! determine different filter bins for the same tally in order to score to it + + TALLY_LOOP: do i = 1, active_cell_to_cell_tallies % size() + ! Get index of tally and pointer to tally + i_tally = active_cell_to_cell_tallies % get_item(i) + t => tallies(i_tally) + + ! Find the first bin in each filter. There may be more than one matching + ! bin per filter, but we'll deal with those later. + do i_filt = 1, size(t % filters) +! print *, "1st bin Matching bins before", matching_bins(i_filt) + call t % filters(i_filt) % obj % get_next_bin(p, t % estimator, & + NO_BIN_FOUND, matching_bins(i_filt), filter_weights(i_filt)) +! print *, "1st Matching bins after", matching_bins(i_filt) + ! If there are no valid bins for this filter, then there is nothing to + ! score and we can move on to the next tally. + if (matching_bins(i_filt) == NO_BIN_FOUND) cycle TALLY_LOOP + end do + + ! ======================================================================== + ! Loop until we've covered all valid bins on each of the filters. + + FILTER_LOOP: do + + ! Determine scoring index and weight for this filter combination + filter_index = sum((matching_bins(1:size(t % filters)) - 1) & + * t % stride) + 1 + filter_weight = product(filter_weights(:size(t % filters))) +! print *, "SCORE_CELL_TO_CELL", filter_index + + ! Determine score for each bin + call score_general(p, t, 0, filter_index, & + 0, dummy_real8, filter_weight) +! call score_general(p, t, (k-1)*t % n_score_bins, filter_index, & +! i_nuclide, atom_density, flux * filter_weight) +! 0 for nuclides, atom_density and flux is set to 1 + + ! ====================================================================== + ! Filter logic + + ! If there are no filters, then we are done. + if (size(t % filters) == 0) exit FILTER_LOOP + + ! Increment the filter bins, starting with the last filter. If we get a + ! NO_BIN_FOUND for the last filter, it means we finished all valid bins + ! for that filter, but next-to-last filter might have more than one + ! valid bin so we need to increment that one as well, and so on. + do i_filt = size(t % filters), 1, -1 +! print *, "2nd bin Matching bins before", matching_bins(i_filt) + + call t % filters(i_filt) % obj % get_next_bin(p, t % estimator, & + matching_bins(i_filt), matching_bins(i_filt), & + filter_weights(i_filt)) +! print *, "2nd bin Matching bins after", matching_bins(i_filt) + + if (matching_bins(i_filt) /= NO_BIN_FOUND) exit + end do + + ! If we got all NO_BIN_FOUNDs, then we have finished all valid bins for + ! each of the filters. Exit the loop. + if (all(matching_bins(:size(t % filters)) == NO_BIN_FOUND)) & + exit FILTER_LOOP + + ! Reset all the filters with NO_BIN_FOUND. This will set them back to + ! their first valid bin. + do i_filt = 1, size(t % filters) + if (matching_bins(i_filt) == NO_BIN_FOUND) then + call t % filters(i_filt) % obj % get_next_bin(p, t % estimator, & + matching_bins(i_filt), matching_bins(i_filt), & + filter_weights(i_filt)) + end if + end do + end do FILTER_LOOP + + ! If the user has specified that we can assume all tallies are spatially + ! separate, this implies that once a tally has been scored to, we needn't + ! check the others. This cuts down on overhead when there are many + ! tallies specified + + if (assume_separate) exit TALLY_LOOP + + end do TALLY_LOOP + + ! Reset tally map positioning + position = 0 + + end subroutine score_cell_to_cell + !=============================================================================== ! SCORE_SURFACE_CURRENT tallies surface crossings in a mesh tally by manually ! determining which mesh surfaces were crossed @@ -4284,7 +4420,10 @@ contains end if elseif (user_tallies(i) % type == TALLY_SURFACE_CURRENT) then call active_current_tallies % add(i_user_tallies + i) + elseif (user_tallies(i) % type == TALLY_CELL_TO_CELL) then + call active_cell_to_cell_tallies % add(i_user_tallies + i) end if + end do end subroutine setup_active_usertallies diff --git a/src/tally_filter.F90 b/src/tally_filter.F90 index efa944d90b..77dff58992 100644 --- a/src/tally_filter.F90 +++ b/src/tally_filter.F90 @@ -70,6 +70,45 @@ module tally_filter procedure :: text_label => text_label_cell procedure :: initialize => initialize_cell end type CellFilter + +!=============================================================================== +! CELLFROMTOFILTER specifies which geometric cells particles exit. +!=============================================================================== + type, extends(TallyFilter) :: CellToCellFilter + integer, allocatable :: cells(:) + type(DictIntInt) :: map + contains + procedure :: get_next_bin => get_next_bin_cell_to_cell + procedure :: to_statepoint => to_statepoint_cell_to_cell + procedure :: text_label => text_label_cell_to_cell + procedure :: initialize => initialize_cell_to_cell + end type CellToCellFilter + +!=============================================================================== +! CELLFROMFILTER specifies which geometric cells particles exit. +!=============================================================================== + type, extends(TallyFilter) :: CellFromFilter + integer, allocatable :: cells(:) + type(DictIntInt) :: map + contains + procedure :: get_next_bin => get_next_bin_cell_from + procedure :: to_statepoint => to_statepoint_cell_from + procedure :: text_label => text_label_cell_from + procedure :: initialize => initialize_cell_from + end type CellFromFilter + +!=============================================================================== +! CELLTOFILTER specifies which geometric cells particles enter. +!=============================================================================== + type, extends(TallyFilter) :: CellToFilter + integer, allocatable :: cells(:) + type(DictIntInt) :: map + contains + procedure :: get_next_bin => get_next_bin_cell_to + procedure :: to_statepoint => to_statepoint_cell_to + procedure :: text_label => text_label_cell_to + procedure :: initialize => initialize_cell_to + end type CellToFilter !=============================================================================== ! DISTRIBCELLFILTER specifies which distributed geometric cells tally events @@ -714,6 +753,282 @@ contains label = "Cell " // to_str(cells(this % cells(bin)) % id) end function text_label_cell +!=============================================================================== +! Cell_to_cell_Filter methods +!=============================================================================== + subroutine get_next_bin_cell_to_cell(this, p, estimator, current_bin, & + next_bin, weight) + + class(CellToCellFilter), intent(in) :: this + type(Particle), intent(in) :: p + integer, intent(in) :: estimator + integer, value, intent(in) :: current_bin + integer, intent(out) :: next_bin + real(8), intent(out) :: weight + + integer :: i, start + + ! Find the coordinate level of the last bin we found. + if (current_bin == NO_BIN_FOUND) then + start = 1 + else + do i = 1, p % n_coord + if (p % coord(i) % cell == this % cells(current_bin)) then + start = i + 1 + exit + end if + end do + end if + + ! Starting one coordinate level deeper, find the next bin. !!!!!!!!!!!!!!!!!! OLD !!!!!!!!!!!!!!!!!!!!!!!!!! + next_bin = NO_BIN_FOUND + weight = ERROR_REAL + do i = start, p % n_coord + if (this % map % has_key(p % coord(i) % cell) .and. & + this % map % has_key(p % last_cell)) then + next_bin = this % map % get_key(p % coord(i) % cell) + weight = ONE + exit + end if + end do + end subroutine get_next_bin_cell_to_cell + + subroutine to_statepoint_cell_to_cell(this, filter_group) + class(CellToCellFilter), intent(in) :: this + integer(HID_T), intent(in) :: filter_group + + integer :: i + integer, allocatable :: cell_ids(:) + + call write_dataset(filter_group, "type", "cell") + call write_dataset(filter_group, "n_bins", this % n_bins) + + allocate(cell_ids(size(this % cells))) + do i = 1, size(this % cells) + cell_ids(i) = cells(this % cells(i)) % id + end do + call write_dataset(filter_group, "bins", cell_ids) + end subroutine to_statepoint_cell_to_cell + + subroutine initialize_cell_to_cell(this) + class(CellToCellFilter), intent(inout) :: this + + integer :: i, id + + ! Convert ids to indices. + do i = 1, this % n_bins + id = this % cells(i) +! print *,this % n_bins, this % cells(i), cell_dict % get_key(id) + if (cell_dict % has_key(id)) then + this % cells(i) = cell_dict % get_key(id) + else + call fatal_error("Could not find cell " // trim(to_str(id)) & + &// " specified on tally filter.") + end if + end do + + ! Generate mapping from cell indices to filter bins. + do i = 1, this % n_bins + call this % map % add_key(this % cells(i), i) + end do + end subroutine initialize_cell_to_cell + + function text_label_cell_to_cell(this, bin) result(label) + class(CellToCellFilter), intent(in) :: this + integer, intent(in) :: bin + character(MAX_LINE_LEN) :: label + + label = "Cell to from" // to_str(cells(this % cells(bin)) % id) + end function text_label_cell_to_cell + +!=============================================================================== +! Cell_From_Filter methods +!=============================================================================== + subroutine get_next_bin_cell_from(this, p, estimator, current_bin, & + next_bin, weight) + + class(CellFromFilter), intent(in) :: this + type(Particle), intent(in) :: p + integer, intent(in) :: estimator + integer, value, intent(in) :: current_bin + integer, intent(out) :: next_bin + real(8), intent(out) :: weight + + integer :: i, start + + ! Find the coordinate level of the last bin we found. +! if (current_bin == NO_BIN_FOUND) then +! start = 1 +! else +! do i = 1, p % last_n_coord !shouldnt be done this way, modify later!! +! ! print *, p % last_cell, this % cells(current_bin) +! if (p % last_cell == this % cells(current_bin)) then +! start = i + 1 +! exit +! end if +! end do +! end if + + ! Starting one coordinate level deeper, find the next bin. + next_bin = NO_BIN_FOUND + weight = ERROR_REAL + weight = ONE + if (current_bin == NO_BIN_FOUND) then +! print *, "Looking in map for cell at", p % last_cell," :",this % map % get_key(p % last_cell) + if (this % map % has_key(p % last_cell)) then +! print *, "Found in map for cell at", p % last_cell," :",this % map % get_key(p % last_cell) + next_bin = this % map % get_key(p % last_cell) + end if + else + next_bin = NO_BIN_FOUND + end if + +! do i = start, p % last_n_coord +! if (this % map % has_key(p % last_cell)) then +! next_bin = this % map % get_key(p % coord(i) % last_cell) !need to use coord(i) to dig one level deeper +! weight = ONE +! exit +! end if +! end do + end subroutine get_next_bin_cell_from + + subroutine to_statepoint_cell_from(this, filter_group) + class(CellFromFilter), intent(in) :: this + integer(HID_T), intent(in) :: filter_group + + integer :: i + integer, allocatable :: cell_ids(:) + + call write_dataset(filter_group, "type", "cell") + call write_dataset(filter_group, "n_bins", this % n_bins) + + allocate(cell_ids(size(this % cells))) + do i = 1, size(this % cells) + cell_ids(i) = cells(this % cells(i)) % id + end do + call write_dataset(filter_group, "bins", cell_ids) + end subroutine to_statepoint_cell_from + + subroutine initialize_cell_from(this) + class(CellFromFilter), intent(inout) :: this + + integer :: i, id + + ! Convert ids to indices. + do i = 1, this % n_bins + id = this % cells(i) +! print *,this % n_bins, this % cells(i), cell_dict % get_key(id) + if (cell_dict % has_key(id)) then + this % cells(i) = cell_dict % get_key(id) + else + call fatal_error("Could not find cell " // trim(to_str(id)) & + &// " specified on tally filter.") + end if + end do + + ! Generate mapping from cell indices to filter bins. + do i = 1, this % n_bins +! print *, "Filling bin/cell map : cell",this % cells(i), "bin", i + call this % map % add_key(this % cells(i), i) + end do + end subroutine initialize_cell_from + + function text_label_cell_from(this, bin) result(label) + class(CellFromFilter), intent(in) :: this + integer, intent(in) :: bin + character(MAX_LINE_LEN) :: label + + label = "Cell to from" // to_str(cells(this % cells(bin)) % id) + end function text_label_cell_from + +!=============================================================================== +! Cell_To_Filter methods +!=============================================================================== + subroutine get_next_bin_cell_to(this, p, estimator, current_bin, & + next_bin, weight) + + class(CellToFilter), intent(in) :: this + type(Particle), intent(in) :: p + integer, intent(in) :: estimator + integer, value, intent(in) :: current_bin + integer, intent(out) :: next_bin + real(8), intent(out) :: weight + + integer :: i, start + + ! Find the coordinate level of the last bin we found. + if (current_bin == NO_BIN_FOUND) then + start = 1 + else + do i = 1, p % n_coord + if (p % coord(i) % cell == this % cells(current_bin)) then + start = i + 1 + exit + end if + end do + end if + + ! Starting one coordinate level deeper, find the next bin. + next_bin = NO_BIN_FOUND + weight = ERROR_REAL + do i = start, p % n_coord + if (this % map % has_key(p % coord(i) % cell)) then + next_bin = this % map % get_key(p % coord(i) % cell) + weight = ONE + exit + end if + end do + end subroutine get_next_bin_cell_to + + subroutine to_statepoint_cell_to(this, filter_group) + class(CellToFilter), intent(in) :: this + integer(HID_T), intent(in) :: filter_group + + integer :: i + integer, allocatable :: cell_ids(:) + + call write_dataset(filter_group, "type", "cell") + call write_dataset(filter_group, "n_bins", this % n_bins) + + allocate(cell_ids(size(this % cells))) + do i = 1, size(this % cells) + cell_ids(i) = cells(this % cells(i)) % id + end do + call write_dataset(filter_group, "bins", cell_ids) + end subroutine to_statepoint_cell_to + + subroutine initialize_cell_to(this) + class(CellToFilter), intent(inout) :: this + + integer :: i, id + + ! Convert ids to indices. + do i = 1, this % n_bins + id = this % cells(i) +! print *,this % n_bins, this % cells(i), cell_dict % get_key(id) + if (cell_dict % has_key(id)) then + this % cells(i) = cell_dict % get_key(id) + else + call fatal_error("Could not find cell " // trim(to_str(id)) & + &// " specified on tally filter.") + end if + end do + + ! Generate mapping from cell indices to filter bins. + do i = 1, this % n_bins + call this % map % add_key(this % cells(i), i) + end do + end subroutine initialize_cell_to + + function text_label_cell_to(this, bin) result(label) + class(CellToFilter), intent(in) :: this + integer, intent(in) :: bin + character(MAX_LINE_LEN) :: label + + label = "Cell to to" // to_str(cells(this % cells(bin)) % id) + end function text_label_cell_to + + !=============================================================================== ! DistribcellFilter methods !=============================================================================== diff --git a/src/tally_initialize.F90 b/src/tally_initialize.F90 index 0e7a50a9c5..4c2c6406c5 100644 --- a/src/tally_initialize.F90 +++ b/src/tally_initialize.F90 @@ -63,6 +63,7 @@ contains t % total_score_bins = t % n_score_bins * t % n_nuclide_bins ! Allocate results array + print *, "Allocating", t % total_score_bins, t % total_filter_bins allocate(t % results(3, t % total_score_bins, t % total_filter_bins)) t % results(:,:,:) = ZERO diff --git a/src/tracking.F90 b/src/tracking.F90 index 7650f598ab..e9fbc86623 100644 --- a/src/tracking.F90 +++ b/src/tracking.F90 @@ -37,6 +37,7 @@ contains integer :: surface_crossed ! surface which particle is on integer :: lattice_translation(3) ! in-lattice translation vector integer :: last_cell ! most recent cell particle was in + integer :: last_n_coord ! most recent number of levels for last_cell integer :: n_event ! number of collisions/crossings real(8) :: d_boundary ! distance to nearest boundary real(8) :: d_collision ! sampled distance to collision @@ -148,7 +149,11 @@ contains ! PARTICLE CROSSES SURFACE if (next_level > 0) p % n_coord = next_level + + ! Saving previous cell data last_cell = p % coord(p % n_coord) % cell + last_n_coord = p % n_coord + p % coord(p % n_coord) % cell = NONE if (any(lattice_translation /= 0)) then ! Particle crosses lattice boundary @@ -158,6 +163,12 @@ contains else ! Particle crosses surface p % surface = surface_crossed + + ! /CHANGE/ Saving stuff on particle + p % last_cell = last_cell + p % last_n_coord = last_n_coord + p % coord(next_level) % last_cell = last_cell + call cross_surface(p, last_cell) p % event = EVENT_SURFACE end if From 21c3d6083cb87d3833ae8bd4023ebd9c27f80555 Mon Sep 17 00:00:00 2001 From: Guillaume Date: Fri, 16 Jun 2017 23:03:27 -0400 Subject: [PATCH 03/54] added projection on normal scores, minor fixes, code clean up --- src/constants.F90 | 5 +- src/endf.F90 | 4 + src/geometry.F90 | 41 +++++---- src/input_xml.F90 | 26 +++--- src/output.F90 | 3 + src/particle_header.F90 | 20 ++--- src/tally.F90 | 79 ++++++++--------- src/tally_filter.F90 | 183 ++++++---------------------------------- src/tracking.F90 | 11 +-- 9 files changed, 120 insertions(+), 252 deletions(-) diff --git a/src/constants.F90 b/src/constants.F90 index 9c52572db6..77d2dcb228 100644 --- a/src/constants.F90 +++ b/src/constants.F90 @@ -305,7 +305,7 @@ module constants EVENT_ABSORB = 2 ! Tally score type - integer, parameter :: N_SCORE_TYPES = 25 + integer, parameter :: N_SCORE_TYPES = 26 integer, parameter :: & SCORE_FLUX = -1, & ! flux SCORE_TOTAL = -2, & ! total reaction rate @@ -331,7 +331,8 @@ module constants 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_CELL_TO_CELL_TYPE = -25 ! cell to cell partial current + SCORE_CELL_TO_CELL = -25, & ! cell to cell partial current + SCORE_CELL_TO_CELL_NORMAL_PROJECTION = -26 ! cell to cell partial current projected on surface normal ! Maximum scattering order supported integer, parameter :: MAX_ANG_ORDER = 10 diff --git a/src/endf.F90 b/src/endf.F90 index 0ba2db388b..f672776a94 100644 --- a/src/endf.F90 +++ b/src/endf.F90 @@ -66,6 +66,10 @@ contains string = "fission-q-prompt" case (SCORE_FISS_Q_RECOV) string = "fission-q-recoverable" + case (SCORE_CELL_TO_CELL) + string = "cell-to-cell-partial-current" + case (SCORE_CELL_TO_CELL_NORMAL_PROJECTION) + string = "cell-to-cell-partial-current-projection-on-normal" ! Normal ENDF-based reactions case (TOTAL_XS) diff --git a/src/geometry.F90 b/src/geometry.F90 index 7283748a0c..03a7de7888 100644 --- a/src/geometry.F90 +++ b/src/geometry.F90 @@ -12,7 +12,7 @@ module geometry use stl_vector, only: VectorInt use string, only: to_str use tally, only: score_surface_current, & - &score_cell_to_cell + &score_partial_current implicit none @@ -399,8 +399,6 @@ contains logical :: rotational ! if rotational periodic BC applied logical :: found ! particle found in universe? class(Surface), pointer :: surf - - !print *, "cross surface called" i_surface = abs(p % surface) surf => surfaces(i_surface)%obj @@ -613,23 +611,33 @@ contains ! cells on the positive side call find_cell(p, found, surf%neighbor_pos) - - ! /CHANGE/ Score cell_to_cell -! print *, "Positive side, calling score_cell_to_cell" - call score_cell_to_cell(p, last_cell) - + + ! Save cosine of angle between surface normal and particle direction + p % normal_proj = dot_product(p % coord(1) % uvw, & + & surf % normal(p % coord(1) % xyz)) / & + & sqrt(dot_product(surf % normal(p % coord(1) % xyz)& + &, surf % normal(p % coord(1) % xyz))) + + ! Score cell to cell partial currents + call score_partial_current(p) + if (found) return elseif (p % surface < 0 .and. allocated(surf%neighbor_neg)) then ! If coming from positive side of surface, search all the neighboring ! cells on the negative side + ! Save cosine of angle between surface normal and particle direction + p % normal_proj = dot_product(p % coord(1) % uvw, & + & surf % normal(p % coord(1) % xyz)) / & + & sqrt(dot_product(surf % normal(p % coord(1) % xyz)& + &, surf % normal(p % coord(1) % xyz))) + call find_cell(p, found, surf%neighbor_neg) - - ! /CHANGE/ Score cell_to_cell -! print *, "Negative side, calling score_cell_to_cell" - call score_cell_to_cell(p, last_cell) - + + ! Score cell to cell partial currents + call score_partial_current(p) + if (found) return end if @@ -662,10 +670,9 @@ contains return end if end if - - ! /CHANGE/ Score cell_to_cell -! print *, "Searched all cells then calling score_cell_to_cell" - call score_cell_to_cell(p, last_cell) + + ! Score cell to cell partial currents + call score_partial_current(p) end subroutine cross_surface diff --git a/src/input_xml.F90 b/src/input_xml.F90 index 7ded1060cb..e74824e1d4 100644 --- a/src/input_xml.F90 +++ b/src/input_xml.F90 @@ -3004,7 +3004,7 @@ contains end if n_words = node_word_count(node_filt, "bins") case ("mesh", "universe", "material", "cell", "distribcell", & - "cellborn", "surface", "delayedgroup") + "cellborn", "cellto", "cellfrom", "surface", "delayedgroup") if (.not. check_for_node(node_filt, "bins")) then call fatal_error("Bins not set in filter " // trim(to_str(filter_id))) end if @@ -3035,7 +3035,7 @@ contains allocate(filt % cells(n_words)) call get_node_array(node_filt, "bins", filt % cells) end select - + case ('cellfrom') ! Allocate and declare the filter type allocate(CellFromFilter :: f % obj) @@ -3046,7 +3046,7 @@ contains allocate(filt % cells(n_words)) call get_node_array(node_filt, "bins", filt % cells) end select - + case ('cellto') ! Allocate and declare the filter type allocate(CellToFilter :: f % obj) @@ -3057,7 +3057,7 @@ contains allocate(filt % cells(n_words)) call get_node_array(node_filt, "bins", filt % cells) end select - + case ('cellborn') ! Allocate and declare the filter type allocate(CellbornFilter :: f % obj) @@ -3414,17 +3414,10 @@ contains t % find_filter(FILTER_DISTRIBCELL) = j type is (CellFilter) t % find_filter(FILTER_CELL) = j - type is (CellFromFilter) t % find_filter(FILTER_CELL_TO_CELL) = j - - t % estimator = ESTIMATOR_ANALOG - type is (CellToFilter) t % find_filter(FILTER_CELL_TO_CELL) = j - - t % estimator = ESTIMATOR_ANALOG - type is (CellbornFilter) t % find_filter(FILTER_CELLBORN) = j type is (MaterialFilter) @@ -3927,13 +3920,14 @@ contains t % find_filter(FILTER_SURFACE) = n_filter t % filter(n_filter) = i_filt - case ('events') - t % score_bins(j) = SCORE_EVENTS - case ('partial_current') t % type = TALLY_CELL_TO_CELL - t % score_bins(j) = SCORE_CELL_TO_CELL_TYPE - + t % score_bins(j) = SCORE_CELL_TO_CELL + case ('projected_partial_current') + t % type = TALLY_CELL_TO_CELL + t % score_bins(j) = SCORE_CELL_TO_CELL_NORMAL_PROJECTION + case ('events') + t % score_bins(j) = SCORE_EVENTS case ('elastic', '(n,elastic)') t % score_bins(j) = ELASTIC case ('(n,2nd)') diff --git a/src/output.F90 b/src/output.F90 index 10d4f0511c..46e8827564 100644 --- a/src/output.F90 +++ b/src/output.F90 @@ -733,6 +733,9 @@ contains score_names(abs(SCORE_INVERSE_VELOCITY)) = "Flux-Weighted Inverse Velocity" score_names(abs(SCORE_FISS_Q_PROMPT)) = "Prompt fission power" score_names(abs(SCORE_FISS_Q_RECOV)) = "Recoverable fission power" + score_names(abs(SCORE_CELL_TO_CELL)) = "Partial current" + score_names(abs(SCORE_CELL_TO_CELL_NORMAL_PROJECTION)) = "Partial & + ¤t projected on normal of surface" ! Create filename for tally output filename = trim(path_output) // "tallies.out" diff --git a/src/particle_header.F90 b/src/particle_header.F90 index c6d478540d..bf22a2c6af 100644 --- a/src/particle_header.F90 +++ b/src/particle_header.F90 @@ -15,9 +15,7 @@ module particle_header !=============================================================================== type LocalCoord - - integer :: last_cell = NONE - + ! Indices in various arrays for this level integer :: cell = NONE integer :: universe = NONE @@ -60,13 +58,11 @@ module particle_header real(8) :: wgt ! particle weight real(8) :: mu ! angle of scatter logical :: alive ! is particle alive? - - ! Previous cell before crossing surface - integer :: last_cell - ! Particle coordinates - integer :: last_n_coord ! number of current coordinates -! type(LocalCoord) :: last_coord(MAX_COORD) ! coordinates for all levels - + + ! Crossing surface data + integer :: last_cell ! last cell the particle was in before crossing + real(8) :: normal_proj ! cos of angle to surface normal when crossing + ! Pre-collision physical data real(8) :: last_xyz_current(3) ! coordinates of the last collision or ! reflective/periodic surface crossing @@ -136,6 +132,7 @@ contains ! clear attributes this % surface = NONE + this % last_cell = NONE this % cell_born = NONE this % material = NONE this % last_material = NONE @@ -143,6 +140,7 @@ contains this % wgt = ONE this % last_wgt = ONE this % absorb_wgt = ZERO + this % normal_proj = 1 this % n_bank = 0 this % wgt_bank = ZERO this % sqrtkT = ERROR_REAL @@ -181,8 +179,6 @@ contains elemental subroutine reset_coord(this) class(LocalCoord), intent(inout) :: this - this % last_cell = NONE - this % cell = NONE this % universe = NONE this % lattice = NONE diff --git a/src/tally.F90 b/src/tally.F90 index 3d18d162f0..ccfd828f15 100644 --- a/src/tally.F90 +++ b/src/tally.F90 @@ -116,12 +116,8 @@ contains !######################################################################### ! Determine appropirate scoring value. - -! print *, "Score bin", score_bin select case(score_bin) - - ! use default case for cell_to_from tallies, with t % estimator == ESTIMATOR_ANALOG case (SCORE_FLUX, SCORE_FLUX_YN) if (t % estimator == ESTIMATOR_ANALOG) then @@ -1180,16 +1176,13 @@ contains end if end if - case(SCORE_CELL_TO_CELL_TYPE) - ! Not using default because score_bin and event_MT thing - ! score_bin is -25, and event_MT seems to match the cell id - ! print*, p % event_MT, score_bin -! print *, "Scoring", p % last_wgt, flux + case(SCORE_CELL_TO_CELL) score = p % last_wgt * flux - + + case(SCORE_CELL_TO_CELL_NORMAL_PROJECTION) + score = p % last_wgt * flux * p % normal_proj + case default - !print *, "Found right score type" - !print * , t % estimator if (t % estimator == ESTIMATOR_ANALOG) then ! Any other score is assumed to be a MT number. Thus, we just need ! to check if it matches the MT number of the event @@ -2091,14 +2084,13 @@ contains case (SCORE_EVENTS) ! Simply count number of scoring events score = ONE - - case(SCORE_CELL_TO_CELL_TYPE) - ! Not using default because score_bin and event_MT thing - ! score_bin is -25, and event_MT seems to match the cell id - ! print*, p % event_MT, score_bin -! print *, "Scoring", p % last_wgt, flux + + case(SCORE_CELL_TO_CELL) score = p % last_wgt * flux + case(SCORE_CELL_TO_CELL_NORMAL_PROJECTION) + score = p % last_wgt * flux * p % normal_proj + end select !######################################################################### @@ -2130,8 +2122,6 @@ contains integer :: num_nm ! Number of N,M orders in harmonic integer :: n ! Moment loop index real(8) :: uvw(3) - -! print *, "Expand and score, score bin", score_bin, "score", score, "filter_index", filter_index select case(score_bin) case (SCORE_SCATTER_N, SCORE_NU_SCATTER_N) @@ -2215,7 +2205,6 @@ contains case default -! print *, "Adding to total in expand-and-score", score !$omp atomic t % results(RESULT_VALUE, score_index, filter_index) = & t % results(RESULT_VALUE, score_index, filter_index) + score @@ -3160,18 +3149,17 @@ contains end subroutine score_collision_tally !=============================================================================== -! SCORE_CELL_TO_CELL tallies partial currents from cell_to to cell_from +! score_partial_current tallies partial currents from cell_to to cell_from !=============================================================================== - subroutine score_cell_to_cell(p, last_cell) + subroutine score_partial_current(p) type(Particle), intent(in) :: p - integer, intent(in) :: last_cell - + ! if particle has leaked - ! DF could preserve leakage aiight - !if p % alive = .false. then - + ! DF could preserve leakage + ! if p % alive = .false. then + !end if integer :: i integer :: i_tally @@ -3183,9 +3171,8 @@ contains real(8) :: flux ! collision estimate of flux real(8) :: filter_weight ! combined weight of all filters logical :: finished ! found all valid bin combinations - real(8) :: dummy_real8 type(TallyObject), pointer :: t - + ! Determine collision estimate of flux if (survival_biasing) then ! We need to account for the fact that some weight was already absorbed @@ -3193,13 +3180,15 @@ contains else flux = p % last_wgt end if - + ! A loop over all tallies is necessary because we need to simultaneously ! determine different filter bins for the same tally in order to score to it + ! However, we can restrict the list to cell_to_cell tallies since we know we + ! are crossing a surface - TALLY_LOOP: do i = 1, active_collision_tallies % size() + TALLY_LOOP: do i = 1, active_cell_to_cell_tallies % size() ! Get index of tally and pointer to tally - i_tally = active_collision_tallies % data(i) + i_tally = active_cell_to_cell_tallies % data(i) t => tallies(i_tally) ! Find all valid bins in each filter if they have not already been found @@ -3245,11 +3234,10 @@ contains filter_weight = filter_weight * filter_matches(i_filt) % weights % & data(i_bin) end do - - + call score_general(p, t, 0, filter_index, & - 0, dummy_real8, filter_weight) - + 0, 0d0, flux * filter_weight) + ! ====================================================================== ! Filter logic @@ -3282,15 +3270,15 @@ contains if (assume_separate) exit TALLY_LOOP end do TALLY_LOOP - + ! Reset filter matches flag filter_matches(:) % bins_present = .false. ! Reset tally map positioning position = 0 - - end subroutine score_cell_to_cell - + + end subroutine score_partial_current + !=============================================================================== ! SCORE_SURFACE_CURRENT tallies surface crossings in a mesh tally by manually ! determining which mesh surfaces were crossed @@ -4584,9 +4572,8 @@ contains call active_current_tallies % push_back(i_user_tallies + i) elseif (user_tallies(i) % type == TALLY_CELL_TO_CELL) then call active_cell_to_cell_tallies % push_back(i_user_tallies + i) - call active_current_tallies % push_back(i_user_tallies + i) end if - + end do call active_tallies % shrink_to_fit() @@ -4594,6 +4581,7 @@ contains call active_tracklength_tallies % shrink_to_fit() call active_collision_tallies % shrink_to_fit() call active_current_tallies % shrink_to_fit() + call active_cell_to_cell_tallies % shrink_to_fit() end subroutine setup_active_usertallies @@ -4617,6 +4605,9 @@ contains else if (active_current_tallies % size() > 0) then call fatal_error("Active current tallies should not exist before CMFD & &tallies!") + else if (active_cell_to_cell_tallies % size() > 0) then + call fatal_error("Active cell to cell tallies should not exist before & + &CMFD tallies!") end if do i = 1, n_cmfd_tallies @@ -4627,6 +4618,7 @@ contains if (cmfd_tallies(i) % type == TALLY_VOLUME) then if (cmfd_tallies(i) % estimator == ESTIMATOR_ANALOG) then call active_analog_tallies % push_back(i_cmfd_tallies + i) + call active_cell_to_cell_tallies % push_back(i_cmfd_tallies + 1) elseif (cmfd_tallies(i) % estimator == ESTIMATOR_TRACKLENGTH) then call active_tracklength_tallies % push_back(i_cmfd_tallies + i) end if @@ -4640,6 +4632,7 @@ contains call active_tracklength_tallies % shrink_to_fit() call active_collision_tallies % shrink_to_fit() call active_current_tallies % shrink_to_fit() + call active_cell_to_cell_tallies % shrink_to_fit() end subroutine setup_active_cmfdtallies diff --git a/src/tally_filter.F90 b/src/tally_filter.F90 index 341e44deb5..b70528bf4e 100644 --- a/src/tally_filter.F90 +++ b/src/tally_filter.F90 @@ -70,20 +70,7 @@ module tally_filter procedure :: text_label => text_label_cell procedure :: initialize => initialize_cell end type CellFilter - -!=============================================================================== -! CELLFROMTOFILTER specifies which geometric cells particles exit. -!=============================================================================== - type, extends(TallyFilter) :: CellToCellFilter - integer, allocatable :: cells(:) - type(DictIntInt) :: map - contains - procedure :: get_next_bin => get_next_bin_cell_to_cell - procedure :: to_statepoint => to_statepoint_cell_to_cell - procedure :: text_label => text_label_cell_to_cell - procedure :: initialize => initialize_cell_to_cell - end type CellToCellFilter - + !=============================================================================== ! CELLFROMFILTER specifies which geometric cells particles exit. !=============================================================================== @@ -96,7 +83,7 @@ module tally_filter procedure :: text_label => text_label_cell_from procedure :: initialize => initialize_cell_from end type CellFromFilter - + !=============================================================================== ! CELLTOFILTER specifies which geometric cells particles enter. !=============================================================================== @@ -757,99 +744,10 @@ contains end function text_label_cell !=============================================================================== -! Cell_to_cell_Filter methods +! CellFromFilter methods !=============================================================================== - subroutine get_next_bin_cell_to_cell(this, p, estimator, current_bin, & + subroutine get_next_bin_cell_from(this, p, estimator, current_bin, & next_bin, weight) - - class(CellToCellFilter), intent(in) :: this - type(Particle), intent(in) :: p - integer, intent(in) :: estimator - integer, value, intent(in) :: current_bin - integer, intent(out) :: next_bin - real(8), intent(out) :: weight - - integer :: i, start - - ! Find the coordinate level of the last bin we found. - if (current_bin == NO_BIN_FOUND) then - start = 1 - else - do i = 1, p % n_coord - if (p % coord(i) % cell == this % cells(current_bin)) then - start = i + 1 - exit - end if - end do - end if - - ! Starting one coordinate level deeper, find the next bin. !!!!!!!!!!!!!!!!!! OLD !!!!!!!!!!!!!!!!!!!!!!!!!! - next_bin = NO_BIN_FOUND - weight = ERROR_REAL - do i = start, p % n_coord - if (this % map % has_key(p % coord(i) % cell) .and. & - this % map % has_key(p % last_cell)) then - next_bin = this % map % get_key(p % coord(i) % cell) - weight = ONE - exit - end if - end do - end subroutine get_next_bin_cell_to_cell - - subroutine to_statepoint_cell_to_cell(this, filter_group) - class(CellToCellFilter), intent(in) :: this - integer(HID_T), intent(in) :: filter_group - - integer :: i - integer, allocatable :: cell_ids(:) - - call write_dataset(filter_group, "type", "cell") - call write_dataset(filter_group, "n_bins", this % n_bins) - - allocate(cell_ids(size(this % cells))) - do i = 1, size(this % cells) - cell_ids(i) = cells(this % cells(i)) % id - end do - call write_dataset(filter_group, "bins", cell_ids) - end subroutine to_statepoint_cell_to_cell - - subroutine initialize_cell_to_cell(this) - class(CellToCellFilter), intent(inout) :: this - - integer :: i, id - - ! Convert ids to indices. - do i = 1, this % n_bins - id = this % cells(i) -! print *,this % n_bins, this % cells(i), cell_dict % get_key(id) - if (cell_dict % has_key(id)) then - this % cells(i) = cell_dict % get_key(id) - else - call fatal_error("Could not find cell " // trim(to_str(id)) & - &// " specified on tally filter.") - end if - end do - - ! Generate mapping from cell indices to filter bins. - do i = 1, this % n_bins - call this % map % add_key(this % cells(i), i) - end do - end subroutine initialize_cell_to_cell - - function text_label_cell_to_cell(this, bin) result(label) - class(CellToCellFilter), intent(in) :: this - integer, intent(in) :: bin - character(MAX_LINE_LEN) :: label - - label = "Cell to from" // to_str(cells(this % cells(bin)) % id) - end function text_label_cell_to_cell - -!=============================================================================== -! Cell_From_Filter methods -!=============================================================================== - subroutine get_next_bin_cell_from(this, p, estimator, current_bin, & - next_bin, weight) - class(CellFromFilter), intent(in) :: this type(Particle), intent(in) :: p integer, intent(in) :: estimator @@ -859,45 +757,24 @@ contains integer :: i, start - ! Find the coordinate level of the last bin we found. -! if (current_bin == NO_BIN_FOUND) then -! start = 1 -! else -! do i = 1, p % last_n_coord !shouldnt be done this way, modify later!! -! ! print *, p % last_cell, this % cells(current_bin) -! if (p % last_cell == this % cells(current_bin)) then -! start = i + 1 -! exit -! end if -! end do -! end if - - ! Starting one coordinate level deeper, find the next bin. + ! Particle can only have one cell_from + ! If current_bin is already a bin, then no need to look for another bin next_bin = NO_BIN_FOUND weight = ERROR_REAL - weight = ONE if (current_bin == NO_BIN_FOUND) then -! print *, "Looking in map for cell at", p % last_cell," :",this % map % get_key(p % last_cell) if (this % map % has_key(p % last_cell)) then -! print *, "Found in map for cell at", p % last_cell," :",this % map % get_key(p % last_cell) next_bin = this % map % get_key(p % last_cell) + weight = ONE end if else next_bin = NO_BIN_FOUND end if - -! do i = start, p % last_n_coord -! if (this % map % has_key(p % last_cell)) then -! next_bin = this % map % get_key(p % coord(i) % last_cell) !need to use coord(i) to dig one level deeper -! weight = ONE -! exit -! end if -! end do + end subroutine get_next_bin_cell_from subroutine to_statepoint_cell_from(this, filter_group) class(CellFromFilter), intent(in) :: this - integer(HID_T), intent(in) :: filter_group + integer(HID_T), intent(in) :: filter_group integer :: i integer, allocatable :: cell_ids(:) @@ -920,7 +797,6 @@ contains ! Convert ids to indices. do i = 1, this % n_bins id = this % cells(i) -! print *,this % n_bins, this % cells(i), cell_dict % get_key(id) if (cell_dict % has_key(id)) then this % cells(i) = cell_dict % get_key(id) else @@ -931,31 +807,29 @@ contains ! Generate mapping from cell indices to filter bins. do i = 1, this % n_bins -! print *, "Filling bin/cell map : cell",this % cells(i), "bin", i call this % map % add_key(this % cells(i), i) end do end subroutine initialize_cell_from - + function text_label_cell_from(this, bin) result(label) class(CellFromFilter), intent(in) :: this - integer, intent(in) :: bin - character(MAX_LINE_LEN) :: label + integer, intent(in) :: bin + character(MAX_LINE_LEN) :: label - label = "Cell to from" // to_str(cells(this % cells(bin)) % id) + label = "Cell from " // to_str(cells(this % cells(bin)) % id) end function text_label_cell_from - + !=============================================================================== -! Cell_To_Filter methods +! CellToFilter methods !=============================================================================== - subroutine get_next_bin_cell_to(this, p, estimator, current_bin, & + subroutine get_next_bin_cell_to(this, p, estimator, current_bin, & next_bin, weight) - class(CellToFilter), intent(in) :: this - type(Particle), intent(in) :: p - integer, intent(in) :: estimator - integer, value, intent(in) :: current_bin - integer, intent(out) :: next_bin - real(8), intent(out) :: weight + type(Particle), intent(in) :: p + integer, intent(in) :: estimator + integer, value, intent(in) :: current_bin + integer, intent(out) :: next_bin + real(8), intent(out) :: weight integer :: i, start @@ -985,7 +859,7 @@ contains subroutine to_statepoint_cell_to(this, filter_group) class(CellToFilter), intent(in) :: this - integer(HID_T), intent(in) :: filter_group + integer(HID_T), intent(in) :: filter_group integer :: i integer, allocatable :: cell_ids(:) @@ -1008,7 +882,7 @@ contains ! Convert ids to indices. do i = 1, this % n_bins id = this % cells(i) -! print *,this % n_bins, this % cells(i), cell_dict % get_key(id) + if (cell_dict % has_key(id)) then this % cells(i) = cell_dict % get_key(id) else @@ -1022,16 +896,15 @@ contains call this % map % add_key(this % cells(i), i) end do end subroutine initialize_cell_to - + function text_label_cell_to(this, bin) result(label) class(CellToFilter), intent(in) :: this - integer, intent(in) :: bin - character(MAX_LINE_LEN) :: label + integer, intent(in) :: bin + character(MAX_LINE_LEN) :: label - label = "Cell to to" // to_str(cells(this % cells(bin)) % id) + label = "Cell to " // to_str(cells(this % cells(bin)) % id) end function text_label_cell_to - - + !=============================================================================== ! DistribcellFilter methods !=============================================================================== diff --git a/src/tracking.F90 b/src/tracking.F90 index e9fbc86623..4a87b4ffa8 100644 --- a/src/tracking.F90 +++ b/src/tracking.F90 @@ -149,11 +149,10 @@ contains ! PARTICLE CROSSES SURFACE if (next_level > 0) p % n_coord = next_level - + ! Saving previous cell data last_cell = p % coord(p % n_coord) % cell - last_n_coord = p % n_coord - + p % coord(p % n_coord) % cell = NONE if (any(lattice_translation /= 0)) then ! Particle crosses lattice boundary @@ -163,12 +162,10 @@ contains else ! Particle crosses surface p % surface = surface_crossed - + ! /CHANGE/ Saving stuff on particle p % last_cell = last_cell - p % last_n_coord = last_n_coord - p % coord(next_level) % last_cell = last_cell - + call cross_surface(p, last_cell) p % event = EVENT_SURFACE end if From 51c698635193aed9a282055d3dd9ded32cb93b7f Mon Sep 17 00:00:00 2001 From: Guillaume Date: Fri, 16 Jun 2017 23:07:59 -0400 Subject: [PATCH 04/54] added filters to python API --- openmc/filter.py | 35 ++++++++++++++++++++++++----------- 1 file changed, 24 insertions(+), 11 deletions(-) diff --git a/openmc/filter.py b/openmc/filter.py index 303687e2f1..449141d77a 100644 --- a/openmc/filter.py +++ b/openmc/filter.py @@ -375,9 +375,7 @@ class Filter(object): """ if filter_bin not in self.bins: - print('Failed to find') - print(filter_bin) - print(self.bins) + msg = 'Unable to get the bin index for Filter since "{0}" ' \ 'is not one of the bins'.format(filter_bin) raise ValueError(msg) @@ -506,7 +504,6 @@ class WithIDFilter(Filter): self._bins = bins - class UniverseFilter(WithIDFilter): """Bins tally event locations based on the Universe they occured in. @@ -600,6 +597,10 @@ class CellFilter(WithIDFilter): @property def bins(self): return self._bins + + @bins.setter + def bins(self, bins): + self._smart_set_bins(bins, openmc.Cell) class CellFromFilter(WithIDFilter): """Bins tally on which couple of cells the neutrons came from. @@ -616,6 +617,8 @@ class CellFromFilter(WithIDFilter): ---------- bins : Integral or Iterable of Integral openmc.Cell IDs. + id : int + Unique identifier for the filter num_bins : Integral The number of filter bins stride : Integral @@ -623,6 +626,13 @@ class CellFromFilter(WithIDFilter): filter's bins. """ + @property + def bins(self): + return self._bins + + @bins.setter + def bins(self, bins): + self._smart_set_bins(bins, openmc.Cell) class CellToFilter(WithIDFilter): """Bins tally on which couple of cells the neutrons went to. @@ -639,6 +649,8 @@ class CellToFilter(WithIDFilter): ---------- bins : Integral or Iterable of Integral openmc.Cell IDs. + id : int + Unique identifier for the filter num_bins : Integral The number of filter bins stride : Integral @@ -646,6 +658,13 @@ class CellToFilter(WithIDFilter): filter's bins. """ + @property + def bins(self): + return self._bins + + @bins.setter + def bins(self, bins): + self._smart_set_bins(bins, openmc.Cell) class CellbornFilter(WithIDFilter): """Bins tally events based on which Cell the neutron was born in. @@ -1066,10 +1085,6 @@ class RealFilter(Filter): def get_bin_index(self, filter_bin): i = np.where(abs(self.bins -filter_bin[1]) < 1e-10)[0] if len(i) == 0: - print('Tally bins', self.bins) - print('Bins searched for:', filter_bin) - print('Search result', i) - print(self.bins - filter_bin[1]) msg = 'Unable to get the bin index for Filter since "{0}" ' \ 'is not one of the bins'.format(filter_bin) raise ValueError(msg) @@ -1111,9 +1126,7 @@ class EnergyFilter(RealFilter): def get_bin_index(self, filter_bin): # Use lower energy bound to find index for RealFilters - - #print(self.bins, filter_bin) - + deltas = np.abs(self.bins - filter_bin[1]) / filter_bin[1] min_delta = np.min(deltas) if min_delta < 1E-3: From fdd5eff2e8ee8572c3e6cbac66f6ce6e3486d0d0 Mon Sep 17 00:00:00 2001 From: Guillaume Date: Fri, 16 Jun 2017 23:28:21 -0400 Subject: [PATCH 05/54] white spaces issues --- src/tally.F90 | 2 +- src/tally_filter.F90 | 8 ++++---- 2 files changed, 5 insertions(+), 5 deletions(-) diff --git a/src/tally.F90 b/src/tally.F90 index ccfd828f15..6cda8d627a 100644 --- a/src/tally.F90 +++ b/src/tally.F90 @@ -2090,7 +2090,7 @@ contains case(SCORE_CELL_TO_CELL_NORMAL_PROJECTION) score = p % last_wgt * flux * p % normal_proj - + end select !######################################################################### diff --git a/src/tally_filter.F90 b/src/tally_filter.F90 index b70528bf4e..194f0422d5 100644 --- a/src/tally_filter.F90 +++ b/src/tally_filter.F90 @@ -746,8 +746,8 @@ contains !=============================================================================== ! CellFromFilter methods !=============================================================================== - subroutine get_next_bin_cell_from(this, p, estimator, current_bin, & - next_bin, weight) + subroutine get_next_bin_cell_from(this, p, estimator, current_bin, & + next_bin, weight) class(CellFromFilter), intent(in) :: this type(Particle), intent(in) :: p integer, intent(in) :: estimator @@ -822,8 +822,8 @@ contains !=============================================================================== ! CellToFilter methods !=============================================================================== - subroutine get_next_bin_cell_to(this, p, estimator, current_bin, & - next_bin, weight) + subroutine get_next_bin_cell_to(this, p, estimator, current_bin, & + next_bin, weight) class(CellToFilter), intent(in) :: this type(Particle), intent(in) :: p integer, intent(in) :: estimator From ab9e589174fc30684c6dd93d355cd466503151ed Mon Sep 17 00:00:00 2001 From: Giud Date: Sun, 18 Jun 2017 15:26:54 -0400 Subject: [PATCH 06/54] deleted some debug prints --- openmc/tallies.py | 23 ++--------------------- 1 file changed, 2 insertions(+), 21 deletions(-) diff --git a/openmc/tallies.py b/openmc/tallies.py index 88e9961927..61dbd4d326 100644 --- a/openmc/tallies.py +++ b/openmc/tallies.py @@ -238,10 +238,7 @@ class Tally(object): def num_filter_bins(self): num_bins = 1 - for self_filter in self.filters: - - #print("Looping through filters, num_bin", self_filter.num_bins, num_bins) - + for self_filter in self.filters: num_bins *= self_filter.num_bins return num_bins @@ -1219,7 +1216,6 @@ class Tally(object): # Find the equivalent Filter in this Tally's list of Filters filter_found = self.find_filter(filter_type) - #print('searched for filter', filter_type, 'and found', filter_found) # Get the index for the requested bin from the Filter and return it filter_index = filter_found.get_bin_index(filter_bin) @@ -1381,12 +1377,6 @@ class Tally(object): else: bins = self_filter.bins - #print('in get_filter_indices') - #print('arguments', filters, filter_bins) - #print('params') - #print(bins) - #print(type(self_filter), self_filter) - # Initialize a NumPy array for the Filter bin indices filter_indices.append(np.zeros(len(bins), dtype=np.int)) @@ -1531,8 +1521,6 @@ class Tally(object): e.g., if the score(s) do not match those in the Tally. """ - #print("\n In get_values") - #print("Tally shape", self.shape) # Ensure that the tally has data if (value == 'mean' and self.mean is None) or \ @@ -1542,13 +1530,6 @@ class Tally(object): (value == 'sum_sq' and self.sum_sq is None): msg = 'The Tally ID="{0}" has no data to return'.format(self.id) raise ValueError(msg) - - #print('\n arguments') - #print(filters) - #print(filter_bins) - #print('\n Tally characteristics') - #print(self) - #print(self.name, self.filters, self.scores,self.num_filter_bins) # Get filter, nuclide and score indices filter_indices = self.get_filter_indices(filters, filter_bins) @@ -1557,7 +1538,7 @@ class Tally(object): # Construct outer product of all three index types with each other indices = np.ix_(filter_indices, nuclide_indices, score_indices) - + # Return the desired result from Tally if value == 'mean': data = self.mean[indices] From 3b46ebdd8915cdd90a019dfcc8d8b31341083a01 Mon Sep 17 00:00:00 2001 From: Giud Date: Sun, 18 Jun 2017 16:56:11 -0400 Subject: [PATCH 07/54] Deleted extra comments --- src/tally.F90 | 9 --------- 1 file changed, 9 deletions(-) diff --git a/src/tally.F90 b/src/tally.F90 index 6cda8d627a..a40e0b80f8 100644 --- a/src/tally.F90 +++ b/src/tally.F90 @@ -3156,10 +3156,6 @@ contains type(Particle), intent(in) :: p - ! if particle has leaked - ! DF could preserve leakage - ! if p % alive = .false. then - !end if integer :: i integer :: i_tally @@ -3181,11 +3177,6 @@ contains flux = p % last_wgt end if - ! A loop over all tallies is necessary because we need to simultaneously - ! determine different filter bins for the same tally in order to score to it - ! However, we can restrict the list to cell_to_cell tallies since we know we - ! are crossing a surface - TALLY_LOOP: do i = 1, active_cell_to_cell_tallies % size() ! Get index of tally and pointer to tally i_tally = active_cell_to_cell_tallies % data(i) From 691173886c0c50da8e00d3949d7dfcdc8bb19e7f Mon Sep 17 00:00:00 2001 From: Giud Date: Sun, 18 Jun 2017 17:37:48 -0400 Subject: [PATCH 08/54] Update tallies.py --- openmc/tallies.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/openmc/tallies.py b/openmc/tallies.py index 61dbd4d326..1641d6164d 100644 --- a/openmc/tallies.py +++ b/openmc/tallies.py @@ -238,7 +238,7 @@ class Tally(object): def num_filter_bins(self): num_bins = 1 - for self_filter in self.filters: + for self_filter in self.filters: num_bins *= self_filter.num_bins return num_bins From 3ab762f617f548039e0ed554271cc3140141835d Mon Sep 17 00:00:00 2001 From: Giud Date: Sun, 18 Jun 2017 17:40:05 -0400 Subject: [PATCH 09/54] Update filter.py --- openmc/filter.py | 5 ++--- 1 file changed, 2 insertions(+), 3 deletions(-) diff --git a/openmc/filter.py b/openmc/filter.py index 449141d77a..6c059d66ee 100644 --- a/openmc/filter.py +++ b/openmc/filter.py @@ -375,7 +375,6 @@ class Filter(object): """ if filter_bin not in self.bins: - msg = 'Unable to get the bin index for Filter since "{0}" ' \ 'is not one of the bins'.format(filter_bin) raise ValueError(msg) @@ -504,6 +503,7 @@ class WithIDFilter(Filter): self._bins = bins + class UniverseFilter(WithIDFilter): """Bins tally event locations based on the Universe they occured in. @@ -1083,7 +1083,7 @@ class RealFilter(Filter): return np.allclose(self.bins, other.bins) def get_bin_index(self, filter_bin): - i = np.where(abs(self.bins -filter_bin[1]) < 1e-10)[0] + i = np.where(i = np.where(self.bins == filter_bin[1])[0] if len(i) == 0: msg = 'Unable to get the bin index for Filter since "{0}" ' \ 'is not one of the bins'.format(filter_bin) @@ -1126,7 +1126,6 @@ class EnergyFilter(RealFilter): def get_bin_index(self, filter_bin): # Use lower energy bound to find index for RealFilters - deltas = np.abs(self.bins - filter_bin[1]) / filter_bin[1] min_delta = np.min(deltas) if min_delta < 1E-3: From 974db1ca5718ae7e3cd155e0d91b0e771e48930d Mon Sep 17 00:00:00 2001 From: Giud Date: Sun, 18 Jun 2017 17:53:05 -0400 Subject: [PATCH 10/54] Update filter.py --- openmc/filter.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/openmc/filter.py b/openmc/filter.py index 6c059d66ee..bc6b0bd35e 100644 --- a/openmc/filter.py +++ b/openmc/filter.py @@ -1083,7 +1083,7 @@ class RealFilter(Filter): return np.allclose(self.bins, other.bins) def get_bin_index(self, filter_bin): - i = np.where(i = np.where(self.bins == filter_bin[1])[0] + i = np.where(self.bins == filter_bin[1])[0] if len(i) == 0: msg = 'Unable to get the bin index for Filter since "{0}" ' \ 'is not one of the bins'.format(filter_bin) From 76ab40293b993e348c8ec0c38baa39f7460a4d96 Mon Sep 17 00:00:00 2001 From: Giud Date: Sun, 18 Jun 2017 17:54:28 -0400 Subject: [PATCH 11/54] deleted extra print --- src/tally_initialize.F90 | 1 - 1 file changed, 1 deletion(-) diff --git a/src/tally_initialize.F90 b/src/tally_initialize.F90 index 2182053c37..46c0283642 100644 --- a/src/tally_initialize.F90 +++ b/src/tally_initialize.F90 @@ -62,7 +62,6 @@ contains t % total_score_bins = t % n_score_bins * t % n_nuclide_bins ! Allocate results array - print *, "Allocating", t % total_score_bins, t % total_filter_bins allocate(t % results(3, t % total_score_bins, t % total_filter_bins)) t % results(:,:,:) = ZERO From 81d122657a4d39b8d371e6af551616cf5e5edf96 Mon Sep 17 00:00:00 2001 From: Giud Date: Sun, 18 Jun 2017 17:56:14 -0400 Subject: [PATCH 12/54] Reworded comment --- src/tracking.F90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/tracking.F90 b/src/tracking.F90 index 4a87b4ffa8..ef0c063717 100644 --- a/src/tracking.F90 +++ b/src/tracking.F90 @@ -163,7 +163,7 @@ contains ! Particle crosses surface p % surface = surface_crossed - ! /CHANGE/ Saving stuff on particle + ! Saving last cell for tallying purposes p % last_cell = last_cell call cross_surface(p, last_cell) From e584ade6848b822f541bde6d89299125ffe73758 Mon Sep 17 00:00:00 2001 From: Guillaume Date: Mon, 19 Jun 2017 20:07:24 -0400 Subject: [PATCH 13/54] implemented half of Sterling s corrections (inheriting cellfilter, moved partial current scoring)+ bug fix on energy filter --- src/constants.F90 | 5 +- src/endf.F90 | 2 - src/geometry.F90 | 23 -------- src/input_xml.F90 | 3 - src/output.F90 | 2 - src/tally.F90 | 35 +++++++----- src/tally_filter.F90 | 128 +++++++++++++++++++++---------------------- src/tracking.F90 | 10 +++- 8 files changed, 95 insertions(+), 113 deletions(-) diff --git a/src/constants.F90 b/src/constants.F90 index 77d2dcb228..4e9fe98654 100644 --- a/src/constants.F90 +++ b/src/constants.F90 @@ -305,7 +305,7 @@ module constants EVENT_ABSORB = 2 ! Tally score type - integer, parameter :: N_SCORE_TYPES = 26 + integer, parameter :: N_SCORE_TYPES = 25 integer, parameter :: & SCORE_FLUX = -1, & ! flux SCORE_TOTAL = -2, & ! total reaction rate @@ -331,8 +331,7 @@ module constants 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_CELL_TO_CELL = -25, & ! cell to cell partial current - SCORE_CELL_TO_CELL_NORMAL_PROJECTION = -26 ! cell to cell partial current projected on surface normal + SCORE_CELL_TO_CELL = -25 ! cell to cell partial current ! Maximum scattering order supported integer, parameter :: MAX_ANG_ORDER = 10 diff --git a/src/endf.F90 b/src/endf.F90 index f672776a94..7f3de44f5f 100644 --- a/src/endf.F90 +++ b/src/endf.F90 @@ -68,8 +68,6 @@ contains string = "fission-q-recoverable" case (SCORE_CELL_TO_CELL) string = "cell-to-cell-partial-current" - case (SCORE_CELL_TO_CELL_NORMAL_PROJECTION) - string = "cell-to-cell-partial-current-projection-on-normal" ! Normal ENDF-based reactions case (TOTAL_XS) diff --git a/src/geometry.F90 b/src/geometry.F90 index 03a7de7888..d6c6bff3b0 100644 --- a/src/geometry.F90 +++ b/src/geometry.F90 @@ -611,33 +611,13 @@ contains ! cells on the positive side call find_cell(p, found, surf%neighbor_pos) - - ! Save cosine of angle between surface normal and particle direction - p % normal_proj = dot_product(p % coord(1) % uvw, & - & surf % normal(p % coord(1) % xyz)) / & - & sqrt(dot_product(surf % normal(p % coord(1) % xyz)& - &, surf % normal(p % coord(1) % xyz))) - - ! Score cell to cell partial currents - call score_partial_current(p) - if (found) return elseif (p % surface < 0 .and. allocated(surf%neighbor_neg)) then ! If coming from positive side of surface, search all the neighboring ! cells on the negative side - ! Save cosine of angle between surface normal and particle direction - p % normal_proj = dot_product(p % coord(1) % uvw, & - & surf % normal(p % coord(1) % xyz)) / & - & sqrt(dot_product(surf % normal(p % coord(1) % xyz)& - &, surf % normal(p % coord(1) % xyz))) - call find_cell(p, found, surf%neighbor_neg) - - ! Score cell to cell partial currents - call score_partial_current(p) - if (found) return end if @@ -671,9 +651,6 @@ contains end if end if - ! Score cell to cell partial currents - call score_partial_current(p) - end subroutine cross_surface !=============================================================================== diff --git a/src/input_xml.F90 b/src/input_xml.F90 index e74824e1d4..ab111dcc5a 100644 --- a/src/input_xml.F90 +++ b/src/input_xml.F90 @@ -3923,9 +3923,6 @@ contains case ('partial_current') t % type = TALLY_CELL_TO_CELL t % score_bins(j) = SCORE_CELL_TO_CELL - case ('projected_partial_current') - t % type = TALLY_CELL_TO_CELL - t % score_bins(j) = SCORE_CELL_TO_CELL_NORMAL_PROJECTION case ('events') t % score_bins(j) = SCORE_EVENTS case ('elastic', '(n,elastic)') diff --git a/src/output.F90 b/src/output.F90 index 46e8827564..bfd5014076 100644 --- a/src/output.F90 +++ b/src/output.F90 @@ -734,8 +734,6 @@ contains score_names(abs(SCORE_FISS_Q_PROMPT)) = "Prompt fission power" score_names(abs(SCORE_FISS_Q_RECOV)) = "Recoverable fission power" score_names(abs(SCORE_CELL_TO_CELL)) = "Partial current" - score_names(abs(SCORE_CELL_TO_CELL_NORMAL_PROJECTION)) = "Partial & - ¤t projected on normal of surface" ! Create filename for tally output filename = trim(path_output) // "tallies.out" diff --git a/src/tally.F90 b/src/tally.F90 index a40e0b80f8..d0f504b4eb 100644 --- a/src/tally.F90 +++ b/src/tally.F90 @@ -119,6 +119,7 @@ contains select case(score_bin) + case (SCORE_FLUX, SCORE_FLUX_YN) if (t % estimator == ESTIMATOR_ANALOG) then ! All events score to a flux bin. We actually use a collision @@ -1176,12 +1177,6 @@ contains end if end if - case(SCORE_CELL_TO_CELL) - score = p % last_wgt * flux - - case(SCORE_CELL_TO_CELL_NORMAL_PROJECTION) - score = p % last_wgt * flux * p % normal_proj - case default if (t % estimator == ESTIMATOR_ANALOG) then ! Any other score is assumed to be a MT number. Thus, we just need @@ -2085,12 +2080,6 @@ contains ! Simply count number of scoring events score = ONE - case(SCORE_CELL_TO_CELL) - score = p % last_wgt * flux - - case(SCORE_CELL_TO_CELL_NORMAL_PROJECTION) - score = p % last_wgt * flux * p % normal_proj - end select !######################################################################### @@ -3161,11 +3150,15 @@ contains integer :: i_tally integer :: i_filt integer :: i_bin + integer :: q ! loop index for scoring bins + integer :: score_bin ! scoring bin, e.g. SCORE_FLUX + integer :: score_index ! scoring bin index integer :: j ! loop index for scoring bins integer :: filter_index ! single index for single bin integer :: matching_bin ! next valid filter bin real(8) :: flux ! collision estimate of flux real(8) :: filter_weight ! combined weight of all filters + real(8) :: score ! analog tally score logical :: finished ! found all valid bin combinations type(TallyObject), pointer :: t @@ -3226,8 +3219,22 @@ contains data(i_bin) end do - call score_general(p, t, 0, filter_index, & - 0, 0d0, flux * filter_weight) + ! Determine score + score = flux * filter_weight + + ! Currently only one score type + SCORE_LOOP: do q = 1, t % n_user_score_bins + + ! determine what type of score bin + score_bin = t % score_bins(q) + + ! determine scoring bin index, no offset from nuclide bins + score_index = q + + ! Expand score if necessary and add to tally results. + call expand_and_score(p, t, score_index, filter_index, score_bin, & + score, q) + end do SCORE_LOOP ! ====================================================================== ! Filter logic diff --git a/src/tally_filter.F90 b/src/tally_filter.F90 index 194f0422d5..eb36e8c9f0 100644 --- a/src/tally_filter.F90 +++ b/src/tally_filter.F90 @@ -87,14 +87,12 @@ module tally_filter !=============================================================================== ! CELLTOFILTER specifies which geometric cells particles enter. !=============================================================================== - type, extends(TallyFilter) :: CellToFilter - integer, allocatable :: cells(:) - type(DictIntInt) :: map + type, extends(CellFilter) :: CellToFilter +! integer, allocatable :: cells(:) +! type(DictIntInt) :: map contains - procedure :: get_next_bin => get_next_bin_cell_to - procedure :: to_statepoint => to_statepoint_cell_to - procedure :: text_label => text_label_cell_to - procedure :: initialize => initialize_cell_to + procedure :: to_statepoint_cell => to_statepoint_cell_to + procedure :: text_label_cell => text_label_cell_to end type CellToFilter !=============================================================================== @@ -779,7 +777,7 @@ contains integer :: i integer, allocatable :: cell_ids(:) - call write_dataset(filter_group, "type", "cell") + call write_dataset(filter_group, "type", "cell from") call write_dataset(filter_group, "n_bins", this % n_bins) allocate(cell_ids(size(this % cells))) @@ -822,40 +820,40 @@ contains !=============================================================================== ! CellToFilter methods !=============================================================================== - subroutine get_next_bin_cell_to(this, p, estimator, current_bin, & - next_bin, weight) - class(CellToFilter), intent(in) :: this - type(Particle), intent(in) :: p - integer, intent(in) :: estimator - integer, value, intent(in) :: current_bin - integer, intent(out) :: next_bin - real(8), intent(out) :: weight - - integer :: i, start - - ! Find the coordinate level of the last bin we found. - if (current_bin == NO_BIN_FOUND) then - start = 1 - else - do i = 1, p % n_coord - if (p % coord(i) % cell == this % cells(current_bin)) then - start = i + 1 - exit - end if - end do - end if - - ! Starting one coordinate level deeper, find the next bin. - next_bin = NO_BIN_FOUND - weight = ERROR_REAL - do i = start, p % n_coord - if (this % map % has_key(p % coord(i) % cell)) then - next_bin = this % map % get_key(p % coord(i) % cell) - weight = ONE - exit - end if - end do - end subroutine get_next_bin_cell_to +! subroutine get_next_bin_cell_to(this, p, estimator, current_bin, & +! next_bin, weight) +! class(CellToFilter), intent(in) :: this +! type(Particle), intent(in) :: p +! integer, intent(in) :: estimator +! integer, value, intent(in) :: current_bin +! integer, intent(out) :: next_bin +! real(8), intent(out) :: weight +! +! integer :: i, start +! +! ! Find the coordinate level of the last bin we found. +! if (current_bin == NO_BIN_FOUND) then +! start = 1 +! else +! do i = 1, p % n_coord +! if (p % coord(i) % cell == this % cells(current_bin)) then +! start = i + 1 +! exit +! end if +! end do +! end if +! +! ! Starting one coordinate level deeper, find the next bin. +! next_bin = NO_BIN_FOUND +! weight = ERROR_REAL +! do i = start, p % n_coord +! if (this % map % has_key(p % coord(i) % cell)) then +! next_bin = this % map % get_key(p % coord(i) % cell) +! weight = ONE +! exit +! end if +! end do +! end subroutine get_next_bin_cell_to subroutine to_statepoint_cell_to(this, filter_group) class(CellToFilter), intent(in) :: this @@ -864,7 +862,7 @@ contains integer :: i integer, allocatable :: cell_ids(:) - call write_dataset(filter_group, "type", "cell") + call write_dataset(filter_group, "type", "cell to") call write_dataset(filter_group, "n_bins", this % n_bins) allocate(cell_ids(size(this % cells))) @@ -874,28 +872,28 @@ contains call write_dataset(filter_group, "bins", cell_ids) end subroutine to_statepoint_cell_to - subroutine initialize_cell_to(this) - class(CellToFilter), intent(inout) :: this - - integer :: i, id - - ! Convert ids to indices. - do i = 1, this % n_bins - id = this % cells(i) - - if (cell_dict % has_key(id)) then - this % cells(i) = cell_dict % get_key(id) - else - call fatal_error("Could not find cell " // trim(to_str(id)) & - &// " specified on tally filter.") - end if - end do - - ! Generate mapping from cell indices to filter bins. - do i = 1, this % n_bins - call this % map % add_key(this % cells(i), i) - end do - end subroutine initialize_cell_to +! subroutine initialize_cell_to(this) +! class(CellToFilter), intent(inout) :: this +! +! integer :: i, id +! +! ! Convert ids to indices. +! do i = 1, this % n_bins +! id = this % cells(i) +! +! if (cell_dict % has_key(id)) then +! this % cells(i) = cell_dict % get_key(id) +! else +! call fatal_error("Could not find cell " // trim(to_str(id)) & +! &// " specified on tally filter.") +! end if +! end do +! +! ! Generate mapping from cell indices to filter bins. +! do i = 1, this % n_bins +! call this % map % add_key(this % cells(i), i) +! end do +! end subroutine initialize_cell_to function text_label_cell_to(this, bin) result(label) class(CellToFilter), intent(in) :: this diff --git a/src/tracking.F90 b/src/tracking.F90 index ef0c063717..d8d0f8a7a5 100644 --- a/src/tracking.F90 +++ b/src/tracking.F90 @@ -15,7 +15,7 @@ module tracking use string, only: to_str use tally, only: score_analog_tally, score_tracklength_tally, & score_collision_tally, score_surface_current, & - score_track_derivative, & + score_track_derivative, score_partial_current, & score_collision_derivative, zero_flux_derivs use track_output, only: initialize_particle_track, write_particle_track, & add_particle_track, finalize_particle_track @@ -166,8 +166,16 @@ contains ! Saving last cell for tallying purposes p % last_cell = last_cell + ! Update last_ information. This is needed to use the same filters as + ! the ones implemented for regular tallies + p % last_uvw = p % coord(p % n_coord) % uvw + p % last_E = p % E + call cross_surface(p, last_cell) p % event = EVENT_SURFACE + + ! Score cell to cell partial currents + call score_partial_current(p) end if else ! ==================================================================== From d952d67fe545ba1d26a1a85fdcfa6206cfbfc9ae Mon Sep 17 00:00:00 2001 From: Guillaume Date: Tue, 20 Jun 2017 08:39:16 -0400 Subject: [PATCH 14/54] made cell_from filter more general --- openmc/filter.py | 3 ++- src/geometry.F90 | 3 +-- src/particle_header.F90 | 55 ++++++++++++++++++++++++++++++++++++----- src/tally_filter.F90 | 49 ++++++++++++++++++++---------------- src/tracking.F90 | 10 +++----- 5 files changed, 83 insertions(+), 37 deletions(-) diff --git a/openmc/filter.py b/openmc/filter.py index bc6b0bd35e..9382c333b3 100644 --- a/openmc/filter.py +++ b/openmc/filter.py @@ -18,7 +18,8 @@ AUTO_FILTER_ID = 10000 _FILTER_TYPES = ['universe', 'material', 'cell', 'cellborn', 'surface', 'mesh', 'energy', 'energyout', 'mu', 'polar', 'azimuthal', - 'distribcell', 'delayedgroup', 'energyfunction', 'celltocell', 'cellfrom', 'cellto'] + 'distribcell', 'delayedgroup', 'energyfunction', 'cellfrom', + 'cellto'] _CURRENT_NAMES = {1: 'x-min out', 2: 'x-min in', 3: 'x-max out', 4: 'x-max in', diff --git a/src/geometry.F90 b/src/geometry.F90 index d6c6bff3b0..f80bd7eb70 100644 --- a/src/geometry.F90 +++ b/src/geometry.F90 @@ -11,8 +11,7 @@ module geometry use surface_header use stl_vector, only: VectorInt use string, only: to_str - use tally, only: score_surface_current, & - &score_partial_current + use tally, only: score_surface_current implicit none diff --git a/src/particle_header.F90 b/src/particle_header.F90 index bf22a2c6af..4f9aa25590 100644 --- a/src/particle_header.F90 +++ b/src/particle_header.F90 @@ -32,6 +32,7 @@ module particle_header logical :: rotated = .false. contains procedure :: reset => reset_coord + procedure :: copy_coord end type LocalCoord !=============================================================================== @@ -48,6 +49,10 @@ module particle_header integer :: n_coord ! number of current coordinates type(LocalCoord) :: coord(MAX_COORD) ! coordinates for all levels + ! Particle coordinates before crossing a surface + integer :: last_n_coord ! number of current coordinates + type(LocalCoord) :: last_coord(MAX_COORD) ! coordinates for all levels + ! Energy Data real(8) :: E ! post-collision energy real(8) :: last_E ! pre-collision energy @@ -59,10 +64,6 @@ module particle_header real(8) :: mu ! angle of scatter logical :: alive ! is particle alive? - ! Crossing surface data - integer :: last_cell ! last cell the particle was in before crossing - real(8) :: normal_proj ! cos of angle to surface normal when crossing - ! Pre-collision physical data real(8) :: last_xyz_current(3) ! coordinates of the last collision or ! reflective/periodic surface crossing @@ -110,6 +111,7 @@ module particle_header procedure :: clear => clear_particle procedure :: initialize_from_source procedure :: create_secondary + procedure :: set_last_coord end type Particle contains @@ -132,7 +134,6 @@ contains ! clear attributes this % surface = NONE - this % last_cell = NONE this % cell_born = NONE this % material = NONE this % last_material = NONE @@ -140,7 +141,6 @@ contains this % wgt = ONE this % last_wgt = ONE this % absorb_wgt = ZERO - this % normal_proj = 1 this % n_bank = 0 this % wgt_bank = ZERO this % sqrtkT = ERROR_REAL @@ -153,6 +153,8 @@ contains ! Set up base level coordinates this % coord(1) % universe = root_universe this % n_coord = 1 + this % last_coord(1) % universe = root_universe + this % last_n_coord = 1 end subroutine initialize_particle @@ -223,6 +225,47 @@ contains end subroutine initialize_from_source +!=============================================================================== +! SET_LAST_COORD copies all coordinate levels from coord to last_coord. This is +! meant to have information about the last cell in the cell_from filter +!=============================================================================== + + subroutine set_last_coord(this) + class(Particle), intent(inout) :: this + integer :: i + + this % last_n_coord = this % n_coord + + ! copy all information at all coordinate levels + do i = 1, this % n_coord + call this % last_coord(i) % copy_coord(this % coord(i)) + end do + + ! reset the rest of former last_coord + do i = this % n_coord + 1, MAX_COORD + call this % last_coord(i) % reset() + end do + + end subroutine set_last_coord + +!=============================================================================== +! COPY_COORD copies one coordinate levels from coord to last_coord. +!=============================================================================== + + elemental subroutine copy_coord(this, coord) + class(LocalCoord), intent(inout) :: this + class(LocalCoord), intent(in) :: coord + + this % cell = coord % cell + this % universe = coord % universe + this % lattice = coord % lattice + this % lattice_x = coord % lattice_x + this % lattice_y = coord % lattice_y + this % lattice_z = coord % lattice_z + this % rotated = coord % rotated + + end subroutine copy_coord + !=============================================================================== ! CREATE_SECONDARY stores the current phase space attributes of the particle in ! the secondary bank and increments the number of sites in the secondary bank. diff --git a/src/tally_filter.F90 b/src/tally_filter.F90 index eb36e8c9f0..c4e72eba5b 100644 --- a/src/tally_filter.F90 +++ b/src/tally_filter.F90 @@ -72,24 +72,21 @@ module tally_filter end type CellFilter !=============================================================================== -! CELLFROMFILTER specifies which geometric cells particles exit. +! CELLFROMFILTER specifies which geometric cells particles exit when crossing a +! surface. !=============================================================================== - type, extends(TallyFilter) :: CellFromFilter - integer, allocatable :: cells(:) - type(DictIntInt) :: map + type, extends(CellFilter) :: CellFromFilter contains - procedure :: get_next_bin => get_next_bin_cell_from - procedure :: to_statepoint => to_statepoint_cell_from - procedure :: text_label => text_label_cell_from - procedure :: initialize => initialize_cell_from + procedure :: get_next_bin_cell => get_next_bin_cell_from + procedure :: to_statepoint_cell => to_statepoint_cell_from + procedure :: text_label_cell => text_label_cell_from end type CellFromFilter !=============================================================================== -! CELLTOFILTER specifies which geometric cells particles enter. +! CELLTOFILTER specifies which geometric cells particles enter when crossing a +! surface. !=============================================================================== type, extends(CellFilter) :: CellToFilter -! integer, allocatable :: cells(:) -! type(DictIntInt) :: map contains procedure :: to_statepoint_cell => to_statepoint_cell_to procedure :: text_label_cell => text_label_cell_to @@ -755,18 +752,28 @@ contains integer :: i, start - ! Particle can only have one cell_from - ! If current_bin is already a bin, then no need to look for another bin + ! Find the coordinate level of the last bin we found. + if (current_bin == NO_BIN_FOUND) then + start = 1 + else + do i = 1, p % last_n_coord + if (p % last_coord(i) % cell == this % cells(current_bin)) then + start = i + 1 + exit + end if + end do + end if + + ! Starting one coordinate level deeper, find the next bin. next_bin = NO_BIN_FOUND weight = ERROR_REAL - if (current_bin == NO_BIN_FOUND) then - if (this % map % has_key(p % last_cell)) then - next_bin = this % map % get_key(p % last_cell) + do i = start, p % last_n_coord + if (this % map % has_key(p % last_coord(i) % cell)) then + next_bin = this % map % get_key(p % last_coord(i) % cell) weight = ONE + exit end if - else - next_bin = NO_BIN_FOUND - end if + end do end subroutine get_next_bin_cell_from @@ -777,7 +784,7 @@ contains integer :: i integer, allocatable :: cell_ids(:) - call write_dataset(filter_group, "type", "cell from") + call write_dataset(filter_group, "type", "cellfrom") call write_dataset(filter_group, "n_bins", this % n_bins) allocate(cell_ids(size(this % cells))) @@ -862,7 +869,7 @@ contains integer :: i integer, allocatable :: cell_ids(:) - call write_dataset(filter_group, "type", "cell to") + call write_dataset(filter_group, "type", "cellto") call write_dataset(filter_group, "n_bins", this % n_bins) allocate(cell_ids(size(this % cells))) diff --git a/src/tracking.F90 b/src/tracking.F90 index d8d0f8a7a5..c05cdf54fb 100644 --- a/src/tracking.F90 +++ b/src/tracking.F90 @@ -37,7 +37,6 @@ contains integer :: surface_crossed ! surface which particle is on integer :: lattice_translation(3) ! in-lattice translation vector integer :: last_cell ! most recent cell particle was in - integer :: last_n_coord ! most recent number of levels for last_cell integer :: n_event ! number of collisions/crossings real(8) :: d_boundary ! distance to nearest boundary real(8) :: d_collision ! sampled distance to collision @@ -152,6 +151,7 @@ contains ! Saving previous cell data last_cell = p % coord(p % n_coord) % cell +! p % set_last_coord() p % coord(p % n_coord) % cell = NONE if (any(lattice_translation /= 0)) then @@ -163,9 +163,6 @@ contains ! Particle crosses surface p % surface = surface_crossed - ! Saving last cell for tallying purposes - p % last_cell = last_cell - ! Update last_ information. This is needed to use the same filters as ! the ones implemented for regular tallies p % last_uvw = p % coord(p % n_coord) % uvw @@ -173,10 +170,9 @@ contains call cross_surface(p, last_cell) p % event = EVENT_SURFACE - - ! Score cell to cell partial currents - call score_partial_current(p) end if + ! Score cell to cell partial currents + if(active_cell_to_cell_tallies%size()>0) call score_partial_current(p) else ! ==================================================================== ! PARTICLE HAS COLLISION From 7edc49bd2b366ff33213bef6f78c950528d6be15 Mon Sep 17 00:00:00 2001 From: Guillaume Date: Tue, 20 Jun 2017 08:41:37 -0400 Subject: [PATCH 15/54] deleted commented code in filter --- src/tally_filter.F90 | 80 -------------------------------------------- 1 file changed, 80 deletions(-) diff --git a/src/tally_filter.F90 b/src/tally_filter.F90 index c4e72eba5b..e19bd1bb40 100644 --- a/src/tally_filter.F90 +++ b/src/tally_filter.F90 @@ -794,28 +794,6 @@ contains call write_dataset(filter_group, "bins", cell_ids) end subroutine to_statepoint_cell_from - subroutine initialize_cell_from(this) - class(CellFromFilter), intent(inout) :: this - - integer :: i, id - - ! Convert ids to indices. - do i = 1, this % n_bins - id = this % cells(i) - if (cell_dict % has_key(id)) then - this % cells(i) = cell_dict % get_key(id) - else - call fatal_error("Could not find cell " // trim(to_str(id)) & - &// " specified on tally filter.") - end if - end do - - ! Generate mapping from cell indices to filter bins. - do i = 1, this % n_bins - call this % map % add_key(this % cells(i), i) - end do - end subroutine initialize_cell_from - function text_label_cell_from(this, bin) result(label) class(CellFromFilter), intent(in) :: this integer, intent(in) :: bin @@ -827,41 +805,6 @@ contains !=============================================================================== ! CellToFilter methods !=============================================================================== -! subroutine get_next_bin_cell_to(this, p, estimator, current_bin, & -! next_bin, weight) -! class(CellToFilter), intent(in) :: this -! type(Particle), intent(in) :: p -! integer, intent(in) :: estimator -! integer, value, intent(in) :: current_bin -! integer, intent(out) :: next_bin -! real(8), intent(out) :: weight -! -! integer :: i, start -! -! ! Find the coordinate level of the last bin we found. -! if (current_bin == NO_BIN_FOUND) then -! start = 1 -! else -! do i = 1, p % n_coord -! if (p % coord(i) % cell == this % cells(current_bin)) then -! start = i + 1 -! exit -! end if -! end do -! end if -! -! ! Starting one coordinate level deeper, find the next bin. -! next_bin = NO_BIN_FOUND -! weight = ERROR_REAL -! do i = start, p % n_coord -! if (this % map % has_key(p % coord(i) % cell)) then -! next_bin = this % map % get_key(p % coord(i) % cell) -! weight = ONE -! exit -! end if -! end do -! end subroutine get_next_bin_cell_to - subroutine to_statepoint_cell_to(this, filter_group) class(CellToFilter), intent(in) :: this integer(HID_T), intent(in) :: filter_group @@ -879,29 +822,6 @@ contains call write_dataset(filter_group, "bins", cell_ids) end subroutine to_statepoint_cell_to -! subroutine initialize_cell_to(this) -! class(CellToFilter), intent(inout) :: this -! -! integer :: i, id -! -! ! Convert ids to indices. -! do i = 1, this % n_bins -! id = this % cells(i) -! -! if (cell_dict % has_key(id)) then -! this % cells(i) = cell_dict % get_key(id) -! else -! call fatal_error("Could not find cell " // trim(to_str(id)) & -! &// " specified on tally filter.") -! end if -! end do -! -! ! Generate mapping from cell indices to filter bins. -! do i = 1, this % n_bins -! call this % map % add_key(this % cells(i), i) -! end do -! end subroutine initialize_cell_to - function text_label_cell_to(this, bin) result(label) class(CellToFilter), intent(in) :: this integer, intent(in) :: bin From a56bbcdbfbbd4ecee3c30a34ee9fefd6a73c9080 Mon Sep 17 00:00:00 2001 From: Guillaume Date: Tue, 20 Jun 2017 08:55:30 -0400 Subject: [PATCH 16/54] white spaces --- src/tally.F90 | 2 +- src/tracking.F90 | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/src/tally.F90 b/src/tally.F90 index d0f504b4eb..80a392c15a 100644 --- a/src/tally.F90 +++ b/src/tally.F90 @@ -3230,7 +3230,7 @@ contains ! determine scoring bin index, no offset from nuclide bins score_index = q - + ! Expand score if necessary and add to tally results. call expand_and_score(p, t, score_index, filter_index, score_bin, & score, q) diff --git a/src/tracking.F90 b/src/tracking.F90 index c05cdf54fb..f70a9e115e 100644 --- a/src/tracking.F90 +++ b/src/tracking.F90 @@ -172,7 +172,7 @@ contains p % event = EVENT_SURFACE end if ! Score cell to cell partial currents - if(active_cell_to_cell_tallies%size()>0) call score_partial_current(p) + if(active_cell_to_cell_tallies%size()>0) call score_partial_current(p) else ! ==================================================================== ! PARTICLE HAS COLLISION From 397d5c4ec8d85bf4b417729d419463fb111a9da0 Mon Sep 17 00:00:00 2001 From: Guillaume Date: Tue, 20 Jun 2017 09:22:43 -0400 Subject: [PATCH 17/54] forgot a call --- src/particle_header.F90 | 1 + src/tracking.F90 | 2 +- 2 files changed, 2 insertions(+), 1 deletion(-) diff --git a/src/particle_header.F90 b/src/particle_header.F90 index 4f9aa25590..903b96bd46 100644 --- a/src/particle_header.F90 +++ b/src/particle_header.F90 @@ -231,6 +231,7 @@ contains !=============================================================================== subroutine set_last_coord(this) + class(Particle), intent(inout) :: this integer :: i diff --git a/src/tracking.F90 b/src/tracking.F90 index f70a9e115e..2bb352acc1 100644 --- a/src/tracking.F90 +++ b/src/tracking.F90 @@ -151,7 +151,7 @@ contains ! Saving previous cell data last_cell = p % coord(p % n_coord) % cell -! p % set_last_coord() + call p % set_last_coord() p % coord(p % n_coord) % cell = NONE if (any(lattice_translation /= 0)) then From 0b20d11a377a3b2d92efe2b8be6c3973246b82aa Mon Sep 17 00:00:00 2001 From: Guillaume Date: Tue, 20 Jun 2017 11:52:02 -0400 Subject: [PATCH 18/54] modified inheriting filter in tally filter --- src/tally_filter.F90 | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/src/tally_filter.F90 b/src/tally_filter.F90 index e19bd1bb40..2ceb31284a 100644 --- a/src/tally_filter.F90 +++ b/src/tally_filter.F90 @@ -77,9 +77,9 @@ module tally_filter !=============================================================================== type, extends(CellFilter) :: CellFromFilter contains - procedure :: get_next_bin_cell => get_next_bin_cell_from - procedure :: to_statepoint_cell => to_statepoint_cell_from - procedure :: text_label_cell => text_label_cell_from + procedure :: get_next_bin => get_next_bin_cell_from + procedure :: to_statepoint => to_statepoint_cell_from + procedure :: text_label => text_label_cell_from end type CellFromFilter !=============================================================================== @@ -88,8 +88,8 @@ module tally_filter !=============================================================================== type, extends(CellFilter) :: CellToFilter contains - procedure :: to_statepoint_cell => to_statepoint_cell_to - procedure :: text_label_cell => text_label_cell_to + procedure :: to_statepoint => to_statepoint_cell_to + procedure :: text_label => text_label_cell_to end type CellToFilter !=============================================================================== From 42d741b46389809f87f293c79db4e1cfaae7791d Mon Sep 17 00:00:00 2001 From: Guillaume Date: Tue, 20 Jun 2017 15:04:04 -0400 Subject: [PATCH 19/54] small fix on whether to look for cell to cell tallies --- src/tracking.F90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/tracking.F90 b/src/tracking.F90 index 2bb352acc1..7935eb3d41 100644 --- a/src/tracking.F90 +++ b/src/tracking.F90 @@ -172,7 +172,7 @@ contains p % event = EVENT_SURFACE end if ! Score cell to cell partial currents - if(active_cell_to_cell_tallies%size()>0) call score_partial_current(p) + if(active_cell_to_cell_tallies%size()>2) call score_partial_current(p) else ! ==================================================================== ! PARTICLE HAS COLLISION From 593cca09b3dc7a086fdf6a46b1277f8f15a9e10b Mon Sep 17 00:00:00 2001 From: Guillaume Date: Tue, 20 Jun 2017 15:04:49 -0400 Subject: [PATCH 20/54] save previous data earlier --- src/tracking.F90 | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/src/tracking.F90 b/src/tracking.F90 index 7935eb3d41..a34b5a28b0 100644 --- a/src/tracking.F90 +++ b/src/tracking.F90 @@ -153,6 +153,11 @@ contains last_cell = p % coord(p % n_coord) % cell call p % set_last_coord() + ! Update last_ information. This is needed to use the same filters as + ! the ones implemented for regular tallies + p % last_uvw = p % coord(p % n_coord) % uvw + p % last_E = p % E + p % coord(p % n_coord) % cell = NONE if (any(lattice_translation /= 0)) then ! Particle crosses lattice boundary @@ -163,11 +168,6 @@ contains ! Particle crosses surface p % surface = surface_crossed - ! Update last_ information. This is needed to use the same filters as - ! the ones implemented for regular tallies - p % last_uvw = p % coord(p % n_coord) % uvw - p % last_E = p % E - call cross_surface(p, last_cell) p % event = EVENT_SURFACE end if From 6ecc40676ee41e8f1b42bfbb1197b66315a5a381 Mon Sep 17 00:00:00 2001 From: Guillaume Date: Tue, 20 Jun 2017 18:06:13 -0400 Subject: [PATCH 21/54] deleted action on cell_to_cell tallies in cmfd tally setup --- src/tally.F90 | 2 -- 1 file changed, 2 deletions(-) diff --git a/src/tally.F90 b/src/tally.F90 index 80a392c15a..3ac7e31657 100644 --- a/src/tally.F90 +++ b/src/tally.F90 @@ -4616,7 +4616,6 @@ contains if (cmfd_tallies(i) % type == TALLY_VOLUME) then if (cmfd_tallies(i) % estimator == ESTIMATOR_ANALOG) then call active_analog_tallies % push_back(i_cmfd_tallies + i) - call active_cell_to_cell_tallies % push_back(i_cmfd_tallies + 1) elseif (cmfd_tallies(i) % estimator == ESTIMATOR_TRACKLENGTH) then call active_tracklength_tallies % push_back(i_cmfd_tallies + i) end if @@ -4630,7 +4629,6 @@ contains call active_tracklength_tallies % shrink_to_fit() call active_collision_tallies % shrink_to_fit() call active_current_tallies % shrink_to_fit() - call active_cell_to_cell_tallies % shrink_to_fit() end subroutine setup_active_cmfdtallies From e61a04ee91cf9ac2c53d1ec0f7adafcc7d3eac56 Mon Sep 17 00:00:00 2001 From: Giud Date: Fri, 23 Jun 2017 19:55:28 -0400 Subject: [PATCH 22/54] Update tracking.F90 --- src/tracking.F90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/tracking.F90 b/src/tracking.F90 index a34b5a28b0..13c05cff08 100644 --- a/src/tracking.F90 +++ b/src/tracking.F90 @@ -172,7 +172,7 @@ contains p % event = EVENT_SURFACE end if ! Score cell to cell partial currents - if(active_cell_to_cell_tallies%size()>2) call score_partial_current(p) + if(active_cell_to_cell_tallies%size()>0) call score_partial_current(p) else ! ==================================================================== ! PARTICLE HAS COLLISION From 345f84b42836ded0a0f9481e204505e643bc2ef2 Mon Sep 17 00:00:00 2001 From: guillaume Date: Mon, 26 Jun 2017 22:01:58 -0400 Subject: [PATCH 23/54] Implemented Paul Romano suggestions, changed score and tally names, used associate, and copied all coord to last%coord --- src/constants.F90 | 9 ++++----- src/endf.F90 | 2 -- src/global.F90 | 4 ++-- src/input_xml.F90 | 19 +++---------------- src/output.F90 | 2 +- src/particle_header.F90 | 25 +------------------------ src/tally.F90 | 27 ++++++++++++++------------- src/tally_filter.F90 | 38 -------------------------------------- src/tracking.F90 | 8 ++++---- 9 files changed, 29 insertions(+), 105 deletions(-) diff --git a/src/constants.F90 b/src/constants.F90 index 4e9fe98654..5684e8c3b9 100644 --- a/src/constants.F90 +++ b/src/constants.F90 @@ -289,7 +289,7 @@ module constants integer, parameter :: & TALLY_VOLUME = 1, & TALLY_SURFACE_CURRENT = 2, & - TALLY_CELL_TO_CELL = 3 + TALLY_SURFACE = 3 ! Tally estimator types integer, parameter :: & @@ -305,7 +305,7 @@ module constants EVENT_ABSORB = 2 ! Tally score type - integer, parameter :: N_SCORE_TYPES = 25 + integer, parameter :: N_SCORE_TYPES = 24 integer, parameter :: & SCORE_FLUX = -1, & ! flux SCORE_TOTAL = -2, & ! total reaction rate @@ -330,8 +330,7 @@ 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_CELL_TO_CELL = -25 ! cell to cell partial current + SCORE_DECAY_RATE = -24 ! delayed neutron precursor decay rate ! Maximum scattering order supported integer, parameter :: MAX_ANG_ORDER = 10 @@ -370,7 +369,7 @@ module constants FILTER_AZIMUTHAL = 12, & FILTER_DELAYEDGROUP = 13, & FILTER_ENERGYFUNCTION = 14, & - FILTER_CELL_TO_CELL = 15 + FILTER_CELLFROM = 15 ! Mesh types integer, parameter :: & diff --git a/src/endf.F90 b/src/endf.F90 index 7f3de44f5f..0ba2db388b 100644 --- a/src/endf.F90 +++ b/src/endf.F90 @@ -66,8 +66,6 @@ contains string = "fission-q-prompt" case (SCORE_FISS_Q_RECOV) string = "fission-q-recoverable" - case (SCORE_CELL_TO_CELL) - string = "cell-to-cell-partial-current" ! Normal ENDF-based reactions case (TOTAL_XS) diff --git a/src/global.F90 b/src/global.F90 index 8fefc11c6a..e484804816 100644 --- a/src/global.F90 +++ b/src/global.F90 @@ -159,7 +159,7 @@ module global type(VectorInt) :: active_current_tallies type(VectorInt) :: active_collision_tallies type(VectorInt) :: active_tallies - type(VectorInt) :: active_cell_to_cell_tallies + type(VectorInt) :: active_surface_tallies ! Global tallies ! 1) collision estimate of k-eff @@ -518,7 +518,7 @@ contains call active_tracklength_tallies % clear() call active_current_tallies % clear() call active_collision_tallies % clear() - call active_cell_to_cell_tallies % clear() + call active_surface_tallies % clear() call active_tallies % clear() ! Deallocate track_identifiers diff --git a/src/input_xml.F90 b/src/input_xml.F90 index ab111dcc5a..cef1bd8655 100644 --- a/src/input_xml.F90 +++ b/src/input_xml.F90 @@ -3047,17 +3047,6 @@ contains call get_node_array(node_filt, "bins", filt % cells) end select - case ('cellto') - ! Allocate and declare the filter type - allocate(CellToFilter :: f % obj) - select type (filt => f % obj) - type is (CellToFilter) - ! Allocate and store bins - filt % n_bins = n_words - allocate(filt % cells(n_words)) - call get_node_array(node_filt, "bins", filt % cells) - end select - case ('cellborn') ! Allocate and declare the filter type allocate(CellbornFilter :: f % obj) @@ -3415,9 +3404,7 @@ contains type is (CellFilter) t % find_filter(FILTER_CELL) = j type is (CellFromFilter) - t % find_filter(FILTER_CELL_TO_CELL) = j - type is (CellToFilter) - t % find_filter(FILTER_CELL_TO_CELL) = j + t % find_filter(FILTER_CELLFROM) = j type is (CellbornFilter) t % find_filter(FILTER_CELLBORN) = j type is (MaterialFilter) @@ -3921,8 +3908,8 @@ contains t % filter(n_filter) = i_filt case ('partial_current') - t % type = TALLY_CELL_TO_CELL - t % score_bins(j) = SCORE_CELL_TO_CELL + t % type = TALLY_SURFACE + t % score_bins(j) = SCORE_CURRENT case ('events') t % score_bins(j) = SCORE_EVENTS case ('elastic', '(n,elastic)') diff --git a/src/output.F90 b/src/output.F90 index bfd5014076..a260cbe032 100644 --- a/src/output.F90 +++ b/src/output.F90 @@ -733,7 +733,7 @@ contains score_names(abs(SCORE_INVERSE_VELOCITY)) = "Flux-Weighted Inverse Velocity" score_names(abs(SCORE_FISS_Q_PROMPT)) = "Prompt fission power" score_names(abs(SCORE_FISS_Q_RECOV)) = "Recoverable fission power" - score_names(abs(SCORE_CELL_TO_CELL)) = "Partial current" + score_names(abs(SCORE_CURRENT)) = "Partial current" ! Create filename for tally output filename = trim(path_output) // "tallies.out" diff --git a/src/particle_header.F90 b/src/particle_header.F90 index 903b96bd46..7fcdf56683 100644 --- a/src/particle_header.F90 +++ b/src/particle_header.F90 @@ -32,7 +32,6 @@ module particle_header logical :: rotated = .false. contains procedure :: reset => reset_coord - procedure :: copy_coord end type LocalCoord !=============================================================================== @@ -236,11 +235,7 @@ contains integer :: i this % last_n_coord = this % n_coord - - ! copy all information at all coordinate levels - do i = 1, this % n_coord - call this % last_coord(i) % copy_coord(this % coord(i)) - end do + this % last_coord = this % coord ! reset the rest of former last_coord do i = this % n_coord + 1, MAX_COORD @@ -249,24 +244,6 @@ contains end subroutine set_last_coord -!=============================================================================== -! COPY_COORD copies one coordinate levels from coord to last_coord. -!=============================================================================== - - elemental subroutine copy_coord(this, coord) - class(LocalCoord), intent(inout) :: this - class(LocalCoord), intent(in) :: coord - - this % cell = coord % cell - this % universe = coord % universe - this % lattice = coord % lattice - this % lattice_x = coord % lattice_x - this % lattice_y = coord % lattice_y - this % lattice_z = coord % lattice_z - this % rotated = coord % rotated - - end subroutine copy_coord - !=============================================================================== ! CREATE_SECONDARY stores the current phase space attributes of the particle in ! the secondary bank and increments the number of sites in the secondary bank. diff --git a/src/tally.F90 b/src/tally.F90 index 3ac7e31657..29bf35b70a 100644 --- a/src/tally.F90 +++ b/src/tally.F90 @@ -3138,14 +3138,14 @@ contains end subroutine score_collision_tally !=============================================================================== -! score_partial_current tallies partial currents from cell_to to cell_from +! score_surface_tally is called at every surface crossing and can be used to +! tally partial currents between two cells !=============================================================================== - subroutine score_partial_current(p) + subroutine score_surface_tally(p) type(Particle), intent(in) :: p - !end if integer :: i integer :: i_tally integer :: i_filt @@ -3160,7 +3160,6 @@ contains real(8) :: filter_weight ! combined weight of all filters real(8) :: score ! analog tally score logical :: finished ! found all valid bin combinations - type(TallyObject), pointer :: t ! Determine collision estimate of flux if (survival_biasing) then @@ -3170,10 +3169,10 @@ contains flux = p % last_wgt end if - TALLY_LOOP: do i = 1, active_cell_to_cell_tallies % size() + TALLY_LOOP: do i = 1, active_surface_tallies % size() ! Get index of tally and pointer to tally - i_tally = active_cell_to_cell_tallies % data(i) - t => tallies(i_tally) + i_tally = active_surface_tallies % data(i) + associate (t => tallies(i_tally)) ! Find all valid bins in each filter if they have not already been found ! for a previous tally. @@ -3233,7 +3232,7 @@ contains ! Expand score if necessary and add to tally results. call expand_and_score(p, t, score_index, filter_index, score_bin, & - score, q) + score, score_index) end do SCORE_LOOP ! ====================================================================== @@ -3260,6 +3259,8 @@ contains end do FILTER_LOOP + end associate + ! If the user has specified that we can assume all tallies are spatially ! separate, this implies that once a tally has been scored to, we needn't ! check the others. This cuts down on overhead when there are many @@ -3275,7 +3276,7 @@ contains ! Reset tally map positioning position = 0 - end subroutine score_partial_current + end subroutine score_surface_tally !=============================================================================== ! SCORE_SURFACE_CURRENT tallies surface crossings in a mesh tally by manually @@ -4568,8 +4569,8 @@ contains end if elseif (user_tallies(i) % type == TALLY_SURFACE_CURRENT) then call active_current_tallies % push_back(i_user_tallies + i) - elseif (user_tallies(i) % type == TALLY_CELL_TO_CELL) then - call active_cell_to_cell_tallies % push_back(i_user_tallies + i) + elseif (user_tallies(i) % type == TALLY_SURFACE) then + call active_surface_tallies % push_back(i_user_tallies + i) end if end do @@ -4579,7 +4580,7 @@ contains call active_tracklength_tallies % shrink_to_fit() call active_collision_tallies % shrink_to_fit() call active_current_tallies % shrink_to_fit() - call active_cell_to_cell_tallies % shrink_to_fit() + call active_surface_tallies % shrink_to_fit() end subroutine setup_active_usertallies @@ -4603,7 +4604,7 @@ contains else if (active_current_tallies % size() > 0) then call fatal_error("Active current tallies should not exist before CMFD & &tallies!") - else if (active_cell_to_cell_tallies % size() > 0) then + else if (active_surface_tallies % size() > 0) then call fatal_error("Active cell to cell tallies should not exist before & &CMFD tallies!") end if diff --git a/src/tally_filter.F90 b/src/tally_filter.F90 index 2ceb31284a..0e6518737a 100644 --- a/src/tally_filter.F90 +++ b/src/tally_filter.F90 @@ -82,16 +82,6 @@ module tally_filter procedure :: text_label => text_label_cell_from end type CellFromFilter -!=============================================================================== -! CELLTOFILTER specifies which geometric cells particles enter when crossing a -! surface. -!=============================================================================== - type, extends(CellFilter) :: CellToFilter - contains - procedure :: to_statepoint => to_statepoint_cell_to - procedure :: text_label => text_label_cell_to - end type CellToFilter - !=============================================================================== ! DISTRIBCELLFILTER specifies which distributed geometric cells tally events ! reside in. @@ -802,34 +792,6 @@ contains label = "Cell from " // to_str(cells(this % cells(bin)) % id) end function text_label_cell_from -!=============================================================================== -! CellToFilter methods -!=============================================================================== - subroutine to_statepoint_cell_to(this, filter_group) - class(CellToFilter), intent(in) :: this - integer(HID_T), intent(in) :: filter_group - - integer :: i - integer, allocatable :: cell_ids(:) - - call write_dataset(filter_group, "type", "cellto") - call write_dataset(filter_group, "n_bins", this % n_bins) - - allocate(cell_ids(size(this % cells))) - do i = 1, size(this % cells) - cell_ids(i) = cells(this % cells(i)) % id - end do - call write_dataset(filter_group, "bins", cell_ids) - end subroutine to_statepoint_cell_to - - function text_label_cell_to(this, bin) result(label) - class(CellToFilter), intent(in) :: this - integer, intent(in) :: bin - character(MAX_LINE_LEN) :: label - - label = "Cell to " // to_str(cells(this % cells(bin)) % id) - end function text_label_cell_to - !=============================================================================== ! DistribcellFilter methods !=============================================================================== diff --git a/src/tracking.F90 b/src/tracking.F90 index 13c05cff08..8e67843e89 100644 --- a/src/tracking.F90 +++ b/src/tracking.F90 @@ -15,7 +15,7 @@ module tracking use string, only: to_str use tally, only: score_analog_tally, score_tracklength_tally, & score_collision_tally, score_surface_current, & - score_track_derivative, score_partial_current, & + score_track_derivative, score_surface_tally, & score_collision_derivative, zero_flux_derivs use track_output, only: initialize_particle_track, write_particle_track, & add_particle_track, finalize_particle_track @@ -153,8 +153,8 @@ contains last_cell = p % coord(p % n_coord) % cell call p % set_last_coord() - ! Update last_ information. This is needed to use the same filters as - ! the ones implemented for regular tallies + ! Update last_ data. This is needed to use the same filters in + ! surface tallies as the ones implemented for regular tallies p % last_uvw = p % coord(p % n_coord) % uvw p % last_E = p % E @@ -172,7 +172,7 @@ contains p % event = EVENT_SURFACE end if ! Score cell to cell partial currents - if(active_cell_to_cell_tallies%size()>0) call score_partial_current(p) + if(active_surface_tallies % size() > 0) call score_surface_tally(p) else ! ==================================================================== ! PARTICLE HAS COLLISION From ca9fbb0c8fe6bbfd1fbe7f05ae8ac37b44d8293c Mon Sep 17 00:00:00 2001 From: guillaume Date: Wed, 28 Jun 2017 18:58:28 -0400 Subject: [PATCH 24/54] changes for documentation in base.rst and tallies.rst --- docs/source/pythonapi/base.rst | 1 + docs/source/usersguide/tallies.rst | 10 ++++++++- openmc/filter.py | 36 +++--------------------------- 3 files changed, 13 insertions(+), 34 deletions(-) diff --git a/docs/source/pythonapi/base.rst b/docs/source/pythonapi/base.rst index dde1c91199..e6cf1e1ccc 100644 --- a/docs/source/pythonapi/base.rst +++ b/docs/source/pythonapi/base.rst @@ -117,6 +117,7 @@ Constructing Tallies openmc.UniverseFilter openmc.MaterialFilter openmc.CellFilter + openmc.CellFromFilter openmc.CellbornFilter openmc.SurfaceFilter openmc.MeshFilter diff --git a/docs/source/usersguide/tallies.rst b/docs/source/usersguide/tallies.rst index 3ee2e99d71..46ad97721b 100644 --- a/docs/source/usersguide/tallies.rst +++ b/docs/source/usersguide/tallies.rst @@ -266,7 +266,15 @@ The following tables show all valid scores: | |particle. Note that this score can only be used if | | |a mesh filter has been specified. Furthermore, it | | |may not be used in conjunction with any other | - | |score. | + | |score. Only energy and mesh filters may be used. | + +----------------------+---------------------------------------------------+ + |partial_current |Partial currents on any surface previously defined | + | |in the geometry. Units are particles per source | + | |particle. Note that this score cannot be used with | + | |a current score or a surface or mesh filter. | + | |It may be used in conjunction with any other score.| + | |Surfaces can be defined with cell_from and cell | + | |(to) filters. | +----------------------+---------------------------------------------------+ |events |Number of scoring events. Units are events per | | |source particle. | diff --git a/openmc/filter.py b/openmc/filter.py index 9382c333b3..e822b9553e 100644 --- a/openmc/filter.py +++ b/openmc/filter.py @@ -603,13 +603,14 @@ class CellFilter(WithIDFilter): def bins(self, bins): self._smart_set_bins(bins, openmc.Cell) + class CellFromFilter(WithIDFilter): - """Bins tally on which couple of cells the neutrons came from. + """Bins tally on which Cell the neutron came from. Parameters ---------- bins : openmc.Cell, Integral, or iterable thereof - The Cells to tally. Either openmc.Cell objects or their + The Cell(s) to tally. Either openmc.Cell objects or their Integral ID numbers can be used. filter_id : int Unique identifier for the filter @@ -634,38 +635,7 @@ class CellFromFilter(WithIDFilter): @bins.setter def bins(self, bins): self._smart_set_bins(bins, openmc.Cell) - -class CellToFilter(WithIDFilter): - """Bins tally on which couple of cells the neutrons went to. - Parameters - ---------- - bins : openmc.Cell, Integral, or iterable thereof - The Cells to tally. Either openmc.Cell objects or their - Integral ID numbers can be used. - filter_id : int - Unique identifier for the filter - - Attributes - ---------- - bins : Integral or Iterable of Integral - openmc.Cell IDs. - id : int - Unique identifier for the filter - num_bins : Integral - The number of filter bins - stride : Integral - The number of filter, nuclide and score bins within each of this - filter's bins. - - """ - @property - def bins(self): - return self._bins - - @bins.setter - def bins(self, bins): - self._smart_set_bins(bins, openmc.Cell) class CellbornFilter(WithIDFilter): """Bins tally events based on which Cell the neutron was born in. From aba768d9fdc960305329d3101011395086b6b97d Mon Sep 17 00:00:00 2001 From: guillaume Date: Wed, 28 Jun 2017 20:42:25 -0400 Subject: [PATCH 25/54] changes to documentation --- docs/source/io_formats/tallies.rst | 20 +++++++++++++++----- src/relaxng/tallies.rnc | 14 ++++++++------ src/relaxng/tallies.rng | 2 ++ 3 files changed, 25 insertions(+), 11 deletions(-) diff --git a/docs/source/io_formats/tallies.rst b/docs/source/io_formats/tallies.rst index 9c0a2515e1..25f40c3704 100644 --- a/docs/source/io_formats/tallies.rst +++ b/docs/source/io_formats/tallies.rst @@ -123,8 +123,8 @@ to the scored values. The ``filter`` element has the following attributes/sub-elements: :type: - The type of the filter. Accepted options are "cell", "cellborn", - "material", "universe", "energy", "energyout", "mu", "polar", + The type of the filter. Accepted options are "cell", "cellfrom", + "cellborn", "material", "universe", "energy", "energyout", "mu", "polar", "azimuthal", "mesh", "distribcell", "delayedgroup", and "energyfunction". @@ -154,14 +154,24 @@ For each filter type, the following table describes what the ``bins`` attribute should be set to: :cell: - A list of unique IDs for cells in which the tally should be accumulated. + A list of cells or unique IDs for cells in which the tally should be + accumulated. + +:cellfrom: + This filter allows the tally to be scored when crossing a surface and the + particle came from a specified cell. A list of cell or cell IDs should be + given. + To tally a partial current from a cell to another, this filter should be + used in combination with a cell filter. + This filter should not be used in combination with a surface or meshfilter. :cellborn: This filter allows the tally to be scored to only when particles were - originally born in a specified cell. A list of cell IDs should be given. + originally born in a specified cell. A list of cell or cell IDs should be + given. :material: - A list of unique IDs for matreials in which the tally should be accumulated. + A list of unique IDs for materials in which the tally should be accumulated. :universe: A list of unique IDs for universes in which the tally should be accumulated. diff --git a/src/relaxng/tallies.rnc b/src/relaxng/tallies.rnc index 4674fd4c03..5ab8971e6c 100644 --- a/src/relaxng/tallies.rnc +++ b/src/relaxng/tallies.rnc @@ -35,12 +35,14 @@ element tallies { element filter { (element id { xsd:int } | attribute id { xsd:int }) & ( - ( (element type { ( "cell" | "cellborn" | "material" | "universe" | - "surface" | "distribcell" | "mesh" | "energy" | "energyout" | "mu" | - "polar" | "azimuthal" | "delayedgroup" | "energyfunction") } | - attribute type { ( "cell" | "cellborn" | "material" | "universe" | - "surface" | "distribcell" | "mesh" | "energy" | "energyout" | "mu" | - "polar" | "azimuthal" | "delayedgroup" | "energyfunction") }) & + ( (element type { ( "cell" | "cellfrom" | "cellborn" | "material" | + "universe" | "surface" | "distribcell" | "mesh" | "energy" | + "energyout" | "mu" | "polar" | "azimuthal" | "delayedgroup" | + "energyfunction") } | + attribute type { ( "cell" | "cellfrom" | "cellborn" | "material" | + "universe" | "surface" | "distribcell" | "mesh" | "energy" | + "energyout" | "mu" | "polar" | "azimuthal" | "delayedgroup" | + "energyfunction") }) & (element bins { list { xsd:double+ } } | attribute bins { list { xsd:double+ } }) ) | diff --git a/src/relaxng/tallies.rng b/src/relaxng/tallies.rng index 1c68e57e83..bad2fdce91 100644 --- a/src/relaxng/tallies.rng +++ b/src/relaxng/tallies.rng @@ -168,6 +168,7 @@ cell + cellfrom cellborn material universe @@ -186,6 +187,7 @@ cell + cellfrom cellborn material universe From 6c247fe28408263cdab8ad33d9cec9ff41a5c0fe Mon Sep 17 00:00:00 2001 From: guillaume Date: Wed, 28 Jun 2017 21:51:50 -0400 Subject: [PATCH 26/54] added test suite for surface tally, with partial currents both ways (fuel to water and water to fuel) and couple of filters relevant to crossing surfaces --- src/tally.F90 | 2 +- tests/test_surface_tally/geometry.xml | 11 ++++++ tests/test_surface_tally/materials.xml | 21 +++++++++++ tests/test_surface_tally/results_true.dat | 36 +++++++++++++++++++ tests/test_surface_tally/settings.xml | 12 +++++++ tests/test_surface_tally/tallies.xml | 34 ++++++++++++++++++ .../test_surface_tally/test_surface_tally.py | 11 ++++++ 7 files changed, 126 insertions(+), 1 deletion(-) create mode 100644 tests/test_surface_tally/geometry.xml create mode 100644 tests/test_surface_tally/materials.xml create mode 100644 tests/test_surface_tally/results_true.dat create mode 100644 tests/test_surface_tally/settings.xml create mode 100644 tests/test_surface_tally/tallies.xml create mode 100644 tests/test_surface_tally/test_surface_tally.py diff --git a/src/tally.F90 b/src/tally.F90 index 29bf35b70a..3e96b91938 100644 --- a/src/tally.F90 +++ b/src/tally.F90 @@ -3138,7 +3138,7 @@ contains end subroutine score_collision_tally !=============================================================================== -! score_surface_tally is called at every surface crossing and can be used to +! score_surface_tally is called at every surface crossing and can be used to ! tally partial currents between two cells !=============================================================================== diff --git a/tests/test_surface_tally/geometry.xml b/tests/test_surface_tally/geometry.xml new file mode 100644 index 0000000000..26f4175c68 --- /dev/null +++ b/tests/test_surface_tally/geometry.xml @@ -0,0 +1,11 @@ + + + + + + + + + + + diff --git a/tests/test_surface_tally/materials.xml b/tests/test_surface_tally/materials.xml new file mode 100644 index 0000000000..506f5f74ae --- /dev/null +++ b/tests/test_surface_tally/materials.xml @@ -0,0 +1,21 @@ + + + + + + + + + + + + + + + + + + + + + diff --git a/tests/test_surface_tally/results_true.dat b/tests/test_surface_tally/results_true.dat new file mode 100644 index 0000000000..2566194529 --- /dev/null +++ b/tests/test_surface_tally/results_true.dat @@ -0,0 +1,36 @@ +k-combined: +1.177844E+00 1.713522E-02 +tally 1: +9.400000E-01 +1.778300E-01 +2.561000E+00 +1.315759E+00 +8.618000E+00 +1.487929E+01 +2.559600E+01 +1.311144E+02 +4.700000E-02 +4.870000E-04 +1.010000E-01 +2.207000E-03 +3.780000E-01 +2.889000E-02 +1.089000E+00 +2.432570E-01 +tally 2: +9.070000E-01 +1.648450E-01 +2.495000E+00 +1.247251E+00 +8.515000E+00 +1.452994E+01 +2.539500E+01 +1.290580E+02 +3.800000E-02 +3.220000E-04 +8.400000E-02 +1.480000E-03 +3.540000E-01 +2.528400E-02 +9.770000E-01 +1.957910E-01 diff --git a/tests/test_surface_tally/settings.xml b/tests/test_surface_tally/settings.xml new file mode 100644 index 0000000000..d97ee454a4 --- /dev/null +++ b/tests/test_surface_tally/settings.xml @@ -0,0 +1,12 @@ + + + eigenvalue + 1000 + 10 + 5 + + + -0.62992 -0.62992 -1 0.62992 0.62992 1 + + + diff --git a/tests/test_surface_tally/tallies.xml b/tests/test_surface_tally/tallies.xml new file mode 100644 index 0000000000..1b69016348 --- /dev/null +++ b/tests/test_surface_tally/tallies.xml @@ -0,0 +1,34 @@ + + + + 10000 + + + 10001 + + + 0.0 4000000.0 20000000.0 + + + 0.0 0.785398163397 3.14159265359 + + + 0.0 0.785398163397 3.14159265359 + + + 10001 + + + 10000 + + + 10003 10004 10000 10001 10002 + partial_current + analog + + + 10005 10006 10000 10001 10002 + partial_current + analog + + diff --git a/tests/test_surface_tally/test_surface_tally.py b/tests/test_surface_tally/test_surface_tally.py new file mode 100644 index 0000000000..b04fcc6eba --- /dev/null +++ b/tests/test_surface_tally/test_surface_tally.py @@ -0,0 +1,11 @@ +#!/usr/bin/env python + +import os +import sys +sys.path.insert(0, os.pardir) +from testing_harness import TestHarness + + +if __name__ == '__main__': + harness = TestHarness('statepoint.10.h5') + harness.main() From 9f0f001a688d9899f88e845a074eed7851b66e58 Mon Sep 17 00:00:00 2001 From: guillaume Date: Thu, 29 Jun 2017 18:03:35 -0400 Subject: [PATCH 27/54] enabling surface tallies in input --- src/input_xml.F90 | 8 -------- 1 file changed, 8 deletions(-) diff --git a/src/input_xml.F90 b/src/input_xml.F90 index cef1bd8655..af759fd4e2 100644 --- a/src/input_xml.F90 +++ b/src/input_xml.F90 @@ -3081,7 +3081,6 @@ contains end select case ('surface') - call fatal_error("Surface filter is not yet supported!") ! Allocate and declare the filter type allocate(SurfaceFilter :: f % obj) select type (filt => f % obj) @@ -3442,13 +3441,6 @@ contains ! Store the filter indices call move_alloc(FROM=temp_filter, TO=t % filter) - ! Check that both cell and surface weren't specified - if (t % find_filter(FILTER_CELL) > 0 .and. & - t % find_filter(FILTER_SURFACE) > 0) then - call fatal_error("Cannot specify both cell and surface filters for & - &tally " // trim(to_str(t % id))) - end if - ! ======================================================================= ! READ DATA FOR NUCLIDES From 096d4551c7df5d8538028751b5db8f3781358d52 Mon Sep 17 00:00:00 2001 From: guillaume Date: Thu, 6 Jul 2017 14:16:36 -0400 Subject: [PATCH 28/54] fixed survival biasing in score_surface_tally --- src/tally.F90 | 9 ++------- 1 file changed, 2 insertions(+), 7 deletions(-) diff --git a/src/tally.F90 b/src/tally.F90 index 3e96b91938..7c3ab70207 100644 --- a/src/tally.F90 +++ b/src/tally.F90 @@ -3161,13 +3161,8 @@ contains real(8) :: score ! analog tally score logical :: finished ! found all valid bin combinations - ! Determine collision estimate of flux - if (survival_biasing) then - ! We need to account for the fact that some weight was already absorbed - flux = (p % last_wgt + p % absorb_wgt) - else - flux = p % last_wgt - end if + ! No collision, so no weight change when survival biasing + flux = p % wgt TALLY_LOOP: do i = 1, active_surface_tallies % size() ! Get index of tally and pointer to tally From db5b05b5d818a56b69fa7177c74cb8ae86e8f089 Mon Sep 17 00:00:00 2001 From: guillaume Date: Thu, 6 Jul 2017 14:19:12 -0400 Subject: [PATCH 29/54] fixed pointer issue in calling expand and score from score_surface_tally --- src/tally.F90 | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/src/tally.F90 b/src/tally.F90 index 7c3ab70207..a5c42e6ace 100644 --- a/src/tally.F90 +++ b/src/tally.F90 @@ -3151,6 +3151,7 @@ contains integer :: i_filt integer :: i_bin integer :: q ! loop index for scoring bins + integer :: k ! working index for expand and score integer :: score_bin ! scoring bin, e.g. SCORE_FLUX integer :: score_index ! scoring bin index integer :: j ! loop index for scoring bins @@ -3217,7 +3218,9 @@ contains score = flux * filter_weight ! Currently only one score type + k = 0 SCORE_LOOP: do q = 1, t % n_user_score_bins + k = k + 1 ! determine what type of score bin score_bin = t % score_bins(q) @@ -3227,7 +3230,7 @@ contains ! Expand score if necessary and add to tally results. call expand_and_score(p, t, score_index, filter_index, score_bin, & - score, score_index) + score, k) end do SCORE_LOOP ! ====================================================================== From 81c2f51839704eaa8c912f659acbc09b5228d35a Mon Sep 17 00:00:00 2001 From: guillaume Date: Thu, 6 Jul 2017 16:37:31 -0400 Subject: [PATCH 30/54] changed from ABS to abs to debug travis? --- src/tally_filter.F90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/tally_filter.F90 b/src/tally_filter.F90 index 0e6518737a..62f5afe0e5 100644 --- a/src/tally_filter.F90 +++ b/src/tally_filter.F90 @@ -962,7 +962,7 @@ contains weight = ERROR_REAL if (current_bin == NO_BIN_FOUND) then do i = 1, this % n_bins - if (p % surface == this % surfaces(i)) then + if (abs(p % surface) == this % surfaces(i)) then next_bin = i weight = ONE exit From f7a85dda546a6813fedd15ccc1771611a9a218ec Mon Sep 17 00:00:00 2001 From: guillaume Date: Thu, 6 Jul 2017 19:38:53 -0400 Subject: [PATCH 31/54] Updated docs with surface tallies --- docs/source/io_formats/tallies.rst | 14 ++++++++++---- docs/source/usersguide/tallies.rst | 4 ++-- 2 files changed, 12 insertions(+), 6 deletions(-) diff --git a/docs/source/io_formats/tallies.rst b/docs/source/io_formats/tallies.rst index 25f40c3704..93c58c3489 100644 --- a/docs/source/io_formats/tallies.rst +++ b/docs/source/io_formats/tallies.rst @@ -124,8 +124,8 @@ attributes/sub-elements: :type: The type of the filter. Accepted options are "cell", "cellfrom", - "cellborn", "material", "universe", "energy", "energyout", "mu", "polar", - "azimuthal", "mesh", "distribcell", "delayedgroup", and + "cellborn", "surface", "material", "universe", "energy", "energyout", "mu", + "polar", "azimuthal", "mesh", "distribcell", "delayedgroup", and "energyfunction". :bins: @@ -157,13 +157,19 @@ should be set to: A list of cells or unique IDs for cells in which the tally should be accumulated. +:surface: + This filter allows the tally to be scored when crossing a surface. A list of + surface IDs should be given.It does not specify in which direction the + surface is crossed, and a cellfrom or a cell filter may be used to tally + partial currents. + :cellfrom: This filter allows the tally to be scored when crossing a surface and the particle came from a specified cell. A list of cell or cell IDs should be given. To tally a partial current from a cell to another, this filter should be - used in combination with a cell filter. - This filter should not be used in combination with a surface or meshfilter. + used in combination with a cell filter, to define the other cell. + This filter should not be used in combination with a meshfilter. :cellborn: This filter allows the tally to be scored to only when particles were diff --git a/docs/source/usersguide/tallies.rst b/docs/source/usersguide/tallies.rst index 46ad97721b..521735a833 100644 --- a/docs/source/usersguide/tallies.rst +++ b/docs/source/usersguide/tallies.rst @@ -271,10 +271,10 @@ The following tables show all valid scores: |partial_current |Partial currents on any surface previously defined | | |in the geometry. Units are particles per source | | |particle. Note that this score cannot be used with | - | |a current score or a surface or mesh filter. | + | |a current score or a mesh filter. | | |It may be used in conjunction with any other score.| | |Surfaces can be defined with cell_from and cell | - | |(to) filters. | + | |(to) filters or surface filters. | +----------------------+---------------------------------------------------+ |events |Number of scoring events. Units are events per | | |source particle. | From 06a99dde5e5a04413c3d53eca171bdf840492ebc Mon Sep 17 00:00:00 2001 From: guillaume Date: Thu, 6 Jul 2017 19:59:58 -0400 Subject: [PATCH 32/54] Changed score name from partial current to total current because of surface_filter behavior --- docs/source/io_formats/tallies.rst | 2 +- docs/source/usersguide/tallies.rst | 2 +- src/constants.F90 | 2 +- src/input_xml.F90 | 2 +- src/tally.F90 | 2 +- src/tally_filter.F90 | 3 +-- 6 files changed, 6 insertions(+), 7 deletions(-) diff --git a/docs/source/io_formats/tallies.rst b/docs/source/io_formats/tallies.rst index 93c58c3489..b0de8b6e6f 100644 --- a/docs/source/io_formats/tallies.rst +++ b/docs/source/io_formats/tallies.rst @@ -159,7 +159,7 @@ should be set to: :surface: This filter allows the tally to be scored when crossing a surface. A list of - surface IDs should be given.It does not specify in which direction the + surface IDs should be given. It does not specify in which direction the surface is crossed, and a cellfrom or a cell filter may be used to tally partial currents. diff --git a/docs/source/usersguide/tallies.rst b/docs/source/usersguide/tallies.rst index 521735a833..f928263842 100644 --- a/docs/source/usersguide/tallies.rst +++ b/docs/source/usersguide/tallies.rst @@ -268,7 +268,7 @@ The following tables show all valid scores: | |may not be used in conjunction with any other | | |score. Only energy and mesh filters may be used. | +----------------------+---------------------------------------------------+ - |partial_current |Partial currents on any surface previously defined | + |total-current |Total currents on any surface previously defined | | |in the geometry. Units are particles per source | | |particle. Note that this score cannot be used with | | |a current score or a mesh filter. | diff --git a/src/constants.F90 b/src/constants.F90 index 5684e8c3b9..b2c34be04a 100644 --- a/src/constants.F90 +++ b/src/constants.F90 @@ -319,7 +319,7 @@ module constants SCORE_FISSION = -10, & ! fission rate SCORE_NU_FISSION = -11, & ! neutron production rate SCORE_KAPPA_FISSION = -12, & ! fission energy production rate - SCORE_CURRENT = -13, & ! partial current + SCORE_CURRENT = -13, & ! current SCORE_FLUX_YN = -14, & ! angular moment of flux SCORE_TOTAL_YN = -15, & ! angular moment of total reaction rate SCORE_SCATTER_YN = -16, & ! angular flux-weighted scattering moment (0:N) diff --git a/src/input_xml.F90 b/src/input_xml.F90 index af759fd4e2..8b8367ca78 100644 --- a/src/input_xml.F90 +++ b/src/input_xml.F90 @@ -3899,7 +3899,7 @@ contains t % find_filter(FILTER_SURFACE) = n_filter t % filter(n_filter) = i_filt - case ('partial_current') + case ('total-current') t % type = TALLY_SURFACE t % score_bins(j) = SCORE_CURRENT case ('events') diff --git a/src/tally.F90 b/src/tally.F90 index a5c42e6ace..2c0f7de067 100644 --- a/src/tally.F90 +++ b/src/tally.F90 @@ -3139,7 +3139,7 @@ contains !=============================================================================== ! score_surface_tally is called at every surface crossing and can be used to -! tally partial currents between two cells +! tally total or partial currents between two cells !=============================================================================== subroutine score_surface_tally(p) diff --git a/src/tally_filter.F90 b/src/tally_filter.F90 index 62f5afe0e5..6dc3310018 100644 --- a/src/tally_filter.F90 +++ b/src/tally_filter.F90 @@ -109,8 +109,7 @@ module tally_filter end type CellbornFilter !=============================================================================== -! SURFACEFILTER is currently not implemented for usual geometric surfaces, but -! it is used as a placeholder for mesh surfaces used in current tallies. +! SURFACEFILTER specifies which surface particles are crossing !=============================================================================== type, extends(TallyFilter) :: SurfaceFilter integer, allocatable :: surfaces(:) From e4c5304f5597f18b8684896cd938cea85ffca55a Mon Sep 17 00:00:00 2001 From: guillaume Date: Thu, 6 Jul 2017 20:07:30 -0400 Subject: [PATCH 33/54] Added how to get partial currents in total current description --- docs/source/usersguide/tallies.rst | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/docs/source/usersguide/tallies.rst b/docs/source/usersguide/tallies.rst index f928263842..35b2ff9ac8 100644 --- a/docs/source/usersguide/tallies.rst +++ b/docs/source/usersguide/tallies.rst @@ -272,9 +272,10 @@ The following tables show all valid scores: | |in the geometry. Units are particles per source | | |particle. Note that this score cannot be used with | | |a current score or a mesh filter. | - | |It may be used in conjunction with any other score.| + | |It may be used along with any other filter. | | |Surfaces can be defined with cell_from and cell | - | |(to) filters or surface filters. | + | |(to) filters or surface filters, partial currents | + | |can be obtained with cell_from and cell filters. | +----------------------+---------------------------------------------------+ |events |Number of scoring events. Units are events per | | |source particle. | From 15bbffc6781b14acaf686939a9ceaf4819382c4e Mon Sep 17 00:00:00 2001 From: guillaume Date: Thu, 6 Jul 2017 20:27:58 -0400 Subject: [PATCH 34/54] Added surface_filter, standalone and with cell_from in test_surface_tally --- tests/test_surface_tally/geometry.xml | 16 +++---- tests/test_surface_tally/results_true.dat | 51 +++++++++++++++++++++++ tests/test_surface_tally/tallies.xml | 30 ++++++++++--- 3 files changed, 83 insertions(+), 14 deletions(-) diff --git a/tests/test_surface_tally/geometry.xml b/tests/test_surface_tally/geometry.xml index 26f4175c68..19eae6333b 100644 --- a/tests/test_surface_tally/geometry.xml +++ b/tests/test_surface_tally/geometry.xml @@ -1,11 +1,11 @@ - - - - - - - - + + + + + + + + diff --git a/tests/test_surface_tally/results_true.dat b/tests/test_surface_tally/results_true.dat index 2566194529..f9f473b19c 100644 --- a/tests/test_surface_tally/results_true.dat +++ b/tests/test_surface_tally/results_true.dat @@ -18,6 +18,23 @@ tally 1: 1.089000E+00 2.432570E-01 tally 2: +9.400000E-01 +1.778300E-01 +2.561000E+00 +1.315759E+00 +8.618000E+00 +1.487929E+01 +2.559600E+01 +1.311144E+02 +4.700000E-02 +4.870000E-04 +1.010000E-01 +2.207000E-03 +3.780000E-01 +2.889000E-02 +1.089000E+00 +2.432570E-01 +tally 3: 9.070000E-01 1.648450E-01 2.495000E+00 @@ -34,3 +51,37 @@ tally 2: 2.528400E-02 9.770000E-01 1.957910E-01 +tally 4: +9.070000E-01 +1.648450E-01 +2.495000E+00 +1.247251E+00 +8.515000E+00 +1.452994E+01 +2.539500E+01 +1.290580E+02 +3.800000E-02 +3.220000E-04 +8.400000E-02 +1.480000E-03 +3.540000E-01 +2.528400E-02 +9.770000E-01 +1.957910E-01 +tally 5: +1.847000E+00 +6.843910E-01 +5.056000E+00 +5.124718E+00 +1.713300E+01 +5.881492E+01 +5.099100E+01 +5.203348E+02 +8.500000E-02 +1.595000E-03 +1.850000E-01 +7.259000E-03 +7.320000E-01 +1.080440E-01 +2.066000E+00 +8.753660E-01 diff --git a/tests/test_surface_tally/tallies.xml b/tests/test_surface_tally/tallies.xml index 1b69016348..fd76a84d22 100644 --- a/tests/test_surface_tally/tallies.xml +++ b/tests/test_surface_tally/tallies.xml @@ -15,20 +15,38 @@ 0.0 0.785398163397 3.14159265359 - + + 1 + + 10001 - + 10000 10003 10004 10000 10001 10002 - partial_current + total-current analog - - 10005 10006 10000 10001 10002 - partial_current + + 10003 10005 10000 10001 10002 + total-current + analog + + + 10006 10007 10000 10001 10002 + total-current + analog + + + 10006 10005 10000 10001 10002 + total-current + analog + + + 10005 10000 10001 10002 + total-current analog From f445f3db6dedff8024ab71654b8d664465722a2a Mon Sep 17 00:00:00 2001 From: guillaume Date: Thu, 13 Jul 2017 22:59:17 -0400 Subject: [PATCH 35/54] correction on documentation --- docs/source/io_formats/tallies.rst | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/docs/source/io_formats/tallies.rst b/docs/source/io_formats/tallies.rst index b0de8b6e6f..62075be933 100644 --- a/docs/source/io_formats/tallies.rst +++ b/docs/source/io_formats/tallies.rst @@ -154,18 +154,18 @@ For each filter type, the following table describes what the ``bins`` attribute should be set to: :cell: - A list of cells or unique IDs for cells in which the tally should be + A list of unique IDs for cells in which the tally should be accumulated. :surface: This filter allows the tally to be scored when crossing a surface. A list of - surface IDs should be given. It does not specify in which direction the - surface is crossed, and a cellfrom or a cell filter may be used to tally - partial currents. + surface IDs should be given. By default, net currents are tallied, and to + tally a partial current from one cell to another, this should be used in + combination with a cell or cell_from filter that defines the other cell. :cellfrom: This filter allows the tally to be scored when crossing a surface and the - particle came from a specified cell. A list of cell or cell IDs should be + particle came from a specified cell. A list of cell IDs should be given. To tally a partial current from a cell to another, this filter should be used in combination with a cell filter, to define the other cell. @@ -173,7 +173,7 @@ should be set to: :cellborn: This filter allows the tally to be scored to only when particles were - originally born in a specified cell. A list of cell or cell IDs should be + originally born in a specified cell. A list of cell IDs should be given. :material: From 53875b703424971c818606229599c11d28fa516f Mon Sep 17 00:00:00 2001 From: guillaume Date: Thu, 13 Jul 2017 23:09:27 -0400 Subject: [PATCH 36/54] made surface tally default to net current, and changed tally_surface_current to tally_mesh_current --- src/constants.F90 | 4 ++-- src/input_xml.F90 | 2 +- src/tally.F90 | 4 ++-- src/tally_filter.F90 | 2 +- 4 files changed, 6 insertions(+), 6 deletions(-) diff --git a/src/constants.F90 b/src/constants.F90 index b2c34be04a..73aca25d29 100644 --- a/src/constants.F90 +++ b/src/constants.F90 @@ -288,8 +288,8 @@ module constants ! Tally type integer, parameter :: & TALLY_VOLUME = 1, & - TALLY_SURFACE_CURRENT = 2, & - TALLY_SURFACE = 3 + TALLY_MESH_CURRENT = 2, & + TALLY_SURFACE = 3 ! Tally estimator types integer, parameter :: & diff --git a/src/input_xml.F90 b/src/input_xml.F90 index 8b8367ca78..2bd6058b1f 100644 --- a/src/input_xml.F90 +++ b/src/input_xml.F90 @@ -3810,7 +3810,7 @@ contains t % score_bins(j) = SCORE_FISS_Q_RECOV case ('current') t % score_bins(j) = SCORE_CURRENT - t % type = TALLY_SURFACE_CURRENT + t % type = TALLY_MESH_CURRENT ! Check to make sure that current is the only desired response ! for this tally diff --git a/src/tally.F90 b/src/tally.F90 index 2c0f7de067..04af19f516 100644 --- a/src/tally.F90 +++ b/src/tally.F90 @@ -4565,7 +4565,7 @@ contains elseif (user_tallies(i) % estimator == ESTIMATOR_COLLISION) then call active_collision_tallies % push_back(i_user_tallies + i) end if - elseif (user_tallies(i) % type == TALLY_SURFACE_CURRENT) then + elseif (user_tallies(i) % type == TALLY_MESH_CURRENT) then call active_current_tallies % push_back(i_user_tallies + i) elseif (user_tallies(i) % type == TALLY_SURFACE) then call active_surface_tallies % push_back(i_user_tallies + i) @@ -4618,7 +4618,7 @@ contains elseif (cmfd_tallies(i) % estimator == ESTIMATOR_TRACKLENGTH) then call active_tracklength_tallies % push_back(i_cmfd_tallies + i) end if - elseif (cmfd_tallies(i) % type == TALLY_SURFACE_CURRENT) then + elseif (cmfd_tallies(i) % type == TALLY_MESH_CURRENT) then call active_current_tallies % push_back(i_cmfd_tallies + i) end if end do diff --git a/src/tally_filter.F90 b/src/tally_filter.F90 index 6dc3310018..0c97d99181 100644 --- a/src/tally_filter.F90 +++ b/src/tally_filter.F90 @@ -963,7 +963,7 @@ contains do i = 1, this % n_bins if (abs(p % surface) == this % surfaces(i)) then next_bin = i - weight = ONE + weight = sign(ONE, p % surface) exit end if end do From de803633fd7e81c17dcd853bea0a6ead402f14ef Mon Sep 17 00:00:00 2001 From: guillaume Date: Thu, 13 Jul 2017 23:16:24 -0400 Subject: [PATCH 37/54] missed a few details on previous commit (function descriptions, some uses of tally_surface_current) --- src/cmfd_input.F90 | 6 +++--- src/output.F90 | 2 +- src/tally_filter.F90 | 2 +- src/trigger.F90 | 10 +++++----- 4 files changed, 10 insertions(+), 10 deletions(-) diff --git a/src/cmfd_input.F90 b/src/cmfd_input.F90 index d57c21beb8..68cd94aa6b 100644 --- a/src/cmfd_input.F90 +++ b/src/cmfd_input.F90 @@ -241,7 +241,7 @@ contains ! There are 3 tally types: ! 1: Only an energy in filter-> flux,total,p1 scatter ! 2: Energy in and energy out filter-> nu-scatter,nu-fission -! 3: Surface current +! 3: Mesh current !=============================================================================== subroutine create_cmfd_tally(root) @@ -425,7 +425,7 @@ contains end select end if - ! Duplicate the mesh filter for the surface current tally since other + ! Duplicate the mesh filter for the mesh current tally since other ! tallies use this filter and we need to change the dimension i_filt = i_filt + 1 allocate(MeshFilter :: filters(i_filt) % obj) @@ -597,7 +597,7 @@ contains ! Set macro bins t % score_bins(1) = SCORE_CURRENT - t % type = TALLY_SURFACE_CURRENT + t % type = TALLY_MESH_CURRENT end if end do diff --git a/src/output.F90 b/src/output.F90 index a260cbe032..c9eecd24cb 100644 --- a/src/output.F90 +++ b/src/output.F90 @@ -792,7 +792,7 @@ contains end if ! Handle surface current tallies separately - if (t % type == TALLY_SURFACE_CURRENT) then + if (t % type == TALLY_MESH_CURRENT) then call write_surface_current(t, unit_tally) cycle end if diff --git a/src/tally_filter.F90 b/src/tally_filter.F90 index 0c97d99181..9068306113 100644 --- a/src/tally_filter.F90 +++ b/src/tally_filter.F90 @@ -963,7 +963,7 @@ contains do i = 1, this % n_bins if (abs(p % surface) == this % surfaces(i)) then next_bin = i - weight = sign(ONE, p % surface) + weight = sign(ONE, real(p % surface, 8)) exit end if end do diff --git a/src/trigger.F90 b/src/trigger.F90 index 9eb39e9cc2..80bdd0e167 100644 --- a/src/trigger.F90 +++ b/src/trigger.F90 @@ -161,8 +161,8 @@ contains trigger % rel_err = ZERO trigger % variance = ZERO - ! Surface current tally triggers require special treatment - if (t % type == TALLY_SURFACE_CURRENT) then + ! Mesh current tally triggers require special treatment + if (t % type == TALLY_MESH_CURRENT) then call compute_tally_current(t, trigger) else @@ -281,7 +281,7 @@ contains !=============================================================================== -! COMPUTE_TALLY_CURRENT computes the current for a surface current tally with +! COMPUTE_TALLY_CURRENT computes the current for a mesh current tally with ! precision trigger(s). !=============================================================================== @@ -301,8 +301,8 @@ contains logical :: print_ebin ! should incoming energy bin be displayed? real(8) :: rel_err = ZERO ! temporary relative error of result real(8) :: std_dev = ZERO ! temporary standard deviration of result - type(TallyObject), pointer :: t ! surface current tally - type(TriggerObject) :: trigger ! surface current tally trigger + type(TallyObject), pointer :: t ! mesh current tally + type(TriggerObject) :: trigger ! mesh current tally trigger type(RegularMesh), pointer :: m ! surface current mesh ! Get pointer to mesh From 9fd8bd3fcf3c417f4dba83557e16139ba1a0d3c9 Mon Sep 17 00:00:00 2001 From: guillaume Date: Thu, 13 Jul 2017 23:19:06 -0400 Subject: [PATCH 38/54] better that way --- src/tally_filter.F90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/tally_filter.F90 b/src/tally_filter.F90 index 9068306113..1e40c409fc 100644 --- a/src/tally_filter.F90 +++ b/src/tally_filter.F90 @@ -963,7 +963,7 @@ contains do i = 1, this % n_bins if (abs(p % surface) == this % surfaces(i)) then next_bin = i - weight = sign(ONE, real(p % surface, 8)) + weight = sign(1, p % surface) exit end if end do From 4d6972f2720064e202075204f102bcf35410cfe9 Mon Sep 17 00:00:00 2001 From: guillaume Date: Fri, 14 Jul 2017 00:17:22 -0400 Subject: [PATCH 39/54] only one score type for meshfilters and surface filters, updated docs, optimized last_cell data stored --- docs/source/io_formats/tallies.rst | 1 + docs/source/usersguide/tallies.rst | 28 ++--- src/geometry.F90 | 5 +- src/input_xml.F90 | 187 ++++++++++++++++------------- src/output.F90 | 2 +- src/particle_header.F90 | 26 +--- src/tally_filter.F90 | 6 +- src/tracking.F90 | 8 +- 8 files changed, 128 insertions(+), 135 deletions(-) diff --git a/docs/source/io_formats/tallies.rst b/docs/source/io_formats/tallies.rst index 62075be933..35a205a353 100644 --- a/docs/source/io_formats/tallies.rst +++ b/docs/source/io_formats/tallies.rst @@ -162,6 +162,7 @@ should be set to: surface IDs should be given. By default, net currents are tallied, and to tally a partial current from one cell to another, this should be used in combination with a cell or cell_from filter that defines the other cell. + This filter should not be used in combination with a meshfilter. :cellfrom: This filter allows the tally to be scored when crossing a surface and the diff --git a/docs/source/usersguide/tallies.rst b/docs/source/usersguide/tallies.rst index 35b2ff9ac8..546dcabac2 100644 --- a/docs/source/usersguide/tallies.rst +++ b/docs/source/usersguide/tallies.rst @@ -261,21 +261,19 @@ The following tables show all valid scores: +----------------------+---------------------------------------------------+ |Score | Description | +======================+===================================================+ - |current |Partial currents on the boundaries of each cell in | - | |a mesh. Units are particles per source | - | |particle. Note that this score can only be used if | - | |a mesh filter has been specified. Furthermore, it | - | |may not be used in conjunction with any other | - | |score. Only energy and mesh filters may be used. | - +----------------------+---------------------------------------------------+ - |total-current |Total currents on any surface previously defined | - | |in the geometry. Units are particles per source | - | |particle. Note that this score cannot be used with | - | |a current score or a mesh filter. | - | |It may be used along with any other filter. | - | |Surfaces can be defined with cell_from and cell | - | |(to) filters or surface filters, partial currents | - | |can be obtained with cell_from and cell filters. | + |current |Used in combination with a mesh filter: | + | |Partial currents on the boundaries of each cell in | + | |a mesh. It may not be used in conjunction with any | + | |other score. Only energy and mesh filters may be | + | |used. | + | |Used in combination with a surface filter: | + | |Net currents on any surface previously defined in | + | |the geometry. It may be used along with any other | + | |filter, except mesh filters. | + | |Surfaces can alternatively be defined with cell | + | |from and cell filters and partial currents can be | + | |obtained by using cell_from and cell filters. | + | |Units are particles per source particle. | +----------------------+---------------------------------------------------+ |events |Number of scoring events. Units are events per | | |source particle. | diff --git a/src/geometry.F90 b/src/geometry.F90 index f80bd7eb70..c451b61278 100644 --- a/src/geometry.F90 +++ b/src/geometry.F90 @@ -384,9 +384,8 @@ contains ! the geometry, is reflected, or crosses into a new lattice or cell !=============================================================================== - subroutine cross_surface(p, last_cell) + subroutine cross_surface(p) type(Particle), intent(inout) :: p - integer, intent(in) :: last_cell ! last cell particle was in real(8) :: u ! x-component of direction real(8) :: v ! y-component of direction @@ -469,7 +468,7 @@ contains p%coord(1)%uvw(:) = [u, v, w] / norm ! Reassign particle's cell and surface - p % coord(1) % cell = last_cell + p % coord(1) % cell = p % last_cell(p % last_n_coord) p % surface = -p % surface ! If a reflective surface is coincident with a lattice or universe diff --git a/src/input_xml.F90 b/src/input_xml.F90 index 2bd6058b1f..e0fa63e0db 100644 --- a/src/input_xml.F90 +++ b/src/input_xml.F90 @@ -3809,99 +3809,114 @@ contains case ('fission-q-recoverable') t % score_bins(j) = SCORE_FISS_Q_RECOV case ('current') - t % score_bins(j) = SCORE_CURRENT - t % type = TALLY_MESH_CURRENT - ! Check to make sure that current is the only desired response - ! for this tally - if (n_words > 1) then - call fatal_error("Cannot tally other scores in the & - &same tally as surface currents") - end if + ! Check which type of current is desired: mesh currents or + ! surface currents + if (t % find_filter(FILTER_SURFACE) > 0 .or. & + &t % find_filter(FILTER_CELL) > 0 .or. & + &t % find_filter(FILTER_CELLFROM) > 0) then - ! Get index of mesh filter - i_filter_mesh = t % filter(t % find_filter(FILTER_MESH)) - - ! Check to make sure mesh filter was specified - if (i_filter_mesh == 0) then - call fatal_error("Cannot tally surface current without a mesh & - &filter.") - end if - - ! Get pointer to mesh - select type(filt => filters(i_filter_mesh) % obj) - type is (MeshFilter) - i_mesh = filt % mesh - m => meshes(i_mesh) - end select - - ! Copy filter indices to temporary array - allocate(temp_filter(size(t % filter) + 1)) - temp_filter(1:size(t % filter)) = t % filter - - ! Move allocation back -- temp_filter becomes deallocated during - ! this call - call move_alloc(FROM=temp_filter, TO=t % filter) - n_filter = size(t % filter) - - ! Extend the filters array so we can add a surface filter and mesh - ! filter - call add_filters(2) - - ! Increment number of user filters - n_user_filters = n_user_filters + 2 - - ! Get index of the new mesh filter - i_filt = n_user_filters - 1 - - ! Duplicate the mesh filter since other tallies might use this - ! filter and we need to change the dimension - allocate(MeshFilter :: filters(i_filt) % obj) - select type(filt => filters(i_filt) % obj) - type is (MeshFilter) - filt % id = i_filt - filt % mesh = i_mesh - - ! We need to increase the dimension by one since we also need - ! currents coming into and out of the boundary mesh cells. - filt % n_bins = product(m % dimension + 1) - - ! Add filter to dictionary - call filter_dict % add_key(filt % id, i_filt) - end select - t % filter(t % find_filter(FILTER_MESH)) = i_filt - - ! Get index of the new surface filter - i_filt = n_user_filters - - ! Add surface filter - allocate(SurfaceFilter :: filters(i_filt) % obj) - select type (filt => filters(i_filt) % obj) - type is (SurfaceFilter) - filt % id = i_filt - filt % n_bins = 4 * m % n_dimension - allocate(filt % surfaces(4 * m % n_dimension)) - if (m % n_dimension == 1) then - filt % surfaces = (/ OUT_LEFT, IN_LEFT, OUT_RIGHT, IN_RIGHT /) - elseif (m % n_dimension == 2) then - filt % surfaces = (/ OUT_LEFT, IN_LEFT, OUT_RIGHT, IN_RIGHT, & - OUT_BACK, IN_BACK, OUT_FRONT, IN_FRONT /) - elseif (m % n_dimension == 3) then - filt % surfaces = (/ OUT_LEFT, IN_LEFT, OUT_RIGHT, IN_RIGHT, & - OUT_BACK, IN_BACK, OUT_FRONT, IN_FRONT, OUT_BOTTOM, & - IN_BOTTOM, OUT_TOP, IN_TOP /) + ! Check to make sure that mesh currents are not desired as well + if (t % find_filter(FILTER_MESH) > 0) then + call fatal_error("Cannot tally other mesh currents & + &in the same tally as surface currents") end if - filt % current = .true. + + t % type = TALLY_SURFACE + t % score_bins(j) = SCORE_CURRENT + + else if (t % find_filter(FILTER_MESH) > 0) then + t % score_bins(j) = SCORE_CURRENT + t % type = TALLY_MESH_CURRENT + + ! Check to make sure that current is the only desired response + ! for this tally + if (n_words > 1) then + call fatal_error("Cannot tally other scores in the & + &same tally as surface currents") + end if + + ! Get index of mesh filter + i_filter_mesh = t % filter(t % find_filter(FILTER_MESH)) + + ! Check to make sure mesh filter was specified + if (i_filter_mesh == 0) then + call fatal_error("Cannot tally surface current without a mesh & + &filter.") + end if + + ! Get pointer to mesh + select type(filt => filters(i_filter_mesh) % obj) + type is (MeshFilter) + i_mesh = filt % mesh + m => meshes(i_mesh) + end select + + ! Copy filter indices to temporary array + allocate(temp_filter(size(t % filter) + 1)) + temp_filter(1:size(t % filter)) = t % filter + + ! Move allocation back -- temp_filter becomes deallocated during + ! this call + call move_alloc(FROM=temp_filter, TO=t % filter) + n_filter = size(t % filter) + + ! Extend the filters array so we can add a surface filter and + ! mesh filter + call add_filters(2) + + ! Increment number of user filters + n_user_filters = n_user_filters + 2 + + ! Get index of the new mesh filter + i_filt = n_user_filters - 1 + + ! Duplicate the mesh filter since other tallies might use this + ! filter and we need to change the dimension + allocate(MeshFilter :: filters(i_filt) % obj) + select type(filt => filters(i_filt) % obj) + type is (MeshFilter) + filt % id = i_filt + filt % mesh = i_mesh + + ! We need to increase the dimension by one since we also need + ! currents coming into and out of the boundary mesh cells. + filt % n_bins = product(m % dimension + 1) ! Add filter to dictionary call filter_dict % add_key(filt % id, i_filt) - end select - t % find_filter(FILTER_SURFACE) = n_filter - t % filter(n_filter) = i_filt + end select + t % filter(t % find_filter(FILTER_MESH)) = i_filt + + ! Get index of the new surface filter + i_filt = n_user_filters + + ! Add surface filter + allocate(SurfaceFilter :: filters(i_filt) % obj) + select type (filt => filters(i_filt) % obj) + type is (SurfaceFilter) + filt % id = i_filt + filt % n_bins = 4 * m % n_dimension + allocate(filt % surfaces(4 * m % n_dimension)) + if (m % n_dimension == 1) then + filt % surfaces = (/ OUT_LEFT, IN_LEFT, OUT_RIGHT, IN_RIGHT /) + elseif (m % n_dimension == 2) then + filt % surfaces = (/ OUT_LEFT, IN_LEFT, OUT_RIGHT, IN_RIGHT, & + OUT_BACK, IN_BACK, OUT_FRONT, IN_FRONT /) + elseif (m % n_dimension == 3) then + filt % surfaces = (/ OUT_LEFT, IN_LEFT, OUT_RIGHT, IN_RIGHT, & + OUT_BACK, IN_BACK, OUT_FRONT, IN_FRONT, OUT_BOTTOM, & + IN_BOTTOM, OUT_TOP, IN_TOP /) + end if + filt % current = .true. + + ! Add filter to dictionary + call filter_dict % add_key(filt % id, i_filt) + end select + t % find_filter(FILTER_SURFACE) = n_filter + t % filter(n_filter) = i_filt + end if - case ('total-current') - t % type = TALLY_SURFACE - t % score_bins(j) = SCORE_CURRENT case ('events') t % score_bins(j) = SCORE_EVENTS case ('elastic', '(n,elastic)') diff --git a/src/output.F90 b/src/output.F90 index c9eecd24cb..ade1a5c05e 100644 --- a/src/output.F90 +++ b/src/output.F90 @@ -733,7 +733,7 @@ contains score_names(abs(SCORE_INVERSE_VELOCITY)) = "Flux-Weighted Inverse Velocity" score_names(abs(SCORE_FISS_Q_PROMPT)) = "Prompt fission power" score_names(abs(SCORE_FISS_Q_RECOV)) = "Recoverable fission power" - score_names(abs(SCORE_CURRENT)) = "Partial current" + score_names(abs(SCORE_CURRENT)) = "Current" ! Create filename for tally output filename = trim(path_output) // "tallies.out" diff --git a/src/particle_header.F90 b/src/particle_header.F90 index 7fcdf56683..76876f021d 100644 --- a/src/particle_header.F90 +++ b/src/particle_header.F90 @@ -49,8 +49,8 @@ module particle_header type(LocalCoord) :: coord(MAX_COORD) ! coordinates for all levels ! Particle coordinates before crossing a surface - integer :: last_n_coord ! number of current coordinates - type(LocalCoord) :: last_coord(MAX_COORD) ! coordinates for all levels + integer :: last_n_coord ! number of current coordinates + integer :: last_cell(MAX_COORD) ! coordinates for all levels ! Energy Data real(8) :: E ! post-collision energy @@ -110,7 +110,6 @@ module particle_header procedure :: clear => clear_particle procedure :: initialize_from_source procedure :: create_secondary - procedure :: set_last_coord end type Particle contains @@ -152,7 +151,6 @@ contains ! Set up base level coordinates this % coord(1) % universe = root_universe this % n_coord = 1 - this % last_coord(1) % universe = root_universe this % last_n_coord = 1 end subroutine initialize_particle @@ -224,26 +222,6 @@ contains end subroutine initialize_from_source -!=============================================================================== -! SET_LAST_COORD copies all coordinate levels from coord to last_coord. This is -! meant to have information about the last cell in the cell_from filter -!=============================================================================== - - subroutine set_last_coord(this) - - class(Particle), intent(inout) :: this - integer :: i - - this % last_n_coord = this % n_coord - this % last_coord = this % coord - - ! reset the rest of former last_coord - do i = this % n_coord + 1, MAX_COORD - call this % last_coord(i) % reset() - end do - - end subroutine set_last_coord - !=============================================================================== ! CREATE_SECONDARY stores the current phase space attributes of the particle in ! the secondary bank and increments the number of sites in the secondary bank. diff --git a/src/tally_filter.F90 b/src/tally_filter.F90 index 1e40c409fc..908b6f1a65 100644 --- a/src/tally_filter.F90 +++ b/src/tally_filter.F90 @@ -746,7 +746,7 @@ contains start = 1 else do i = 1, p % last_n_coord - if (p % last_coord(i) % cell == this % cells(current_bin)) then + if (p % last_cell(i) == this % cells(current_bin)) then start = i + 1 exit end if @@ -757,8 +757,8 @@ contains next_bin = NO_BIN_FOUND weight = ERROR_REAL do i = start, p % last_n_coord - if (this % map % has_key(p % last_coord(i) % cell)) then - next_bin = this % map % get_key(p % last_coord(i) % cell) + if (this % map % has_key(p % last_cell(i))) then + next_bin = this % map % get_key(p % last_cell(i)) weight = ONE exit end if diff --git a/src/tracking.F90 b/src/tracking.F90 index 8e67843e89..042e63c27f 100644 --- a/src/tracking.F90 +++ b/src/tracking.F90 @@ -150,8 +150,10 @@ contains if (next_level > 0) p % n_coord = next_level ! Saving previous cell data - last_cell = p % coord(p % n_coord) % cell - call p % set_last_coord() + do j = 1, p % n_coord + p % last_cell(j) = p % coord(j) % cell + end do + p % last_n_coord = p % n_coord ! Update last_ data. This is needed to use the same filters in ! surface tallies as the ones implemented for regular tallies @@ -168,7 +170,7 @@ contains ! Particle crosses surface p % surface = surface_crossed - call cross_surface(p, last_cell) + call cross_surface(p) p % event = EVENT_SURFACE end if ! Score cell to cell partial currents From 68ad95ac95b18319eccd434e989161a06381b40d Mon Sep 17 00:00:00 2001 From: guillaume Date: Fri, 14 Jul 2017 14:22:16 -0400 Subject: [PATCH 40/54] doc change --- docs/source/usersguide/tallies.rst | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/docs/source/usersguide/tallies.rst b/docs/source/usersguide/tallies.rst index 546dcabac2..3a742b59e8 100644 --- a/docs/source/usersguide/tallies.rst +++ b/docs/source/usersguide/tallies.rst @@ -271,8 +271,8 @@ The following tables show all valid scores: | |the geometry. It may be used along with any other | | |filter, except mesh filters. | | |Surfaces can alternatively be defined with cell | - | |from and cell filters and partial currents can be | - | |obtained by using cell_from and cell filters. | + | |from and cell filters thereby resulting in tallying| + | |partial currents. | | |Units are particles per source particle. | +----------------------+---------------------------------------------------+ |events |Number of scoring events. Units are events per | From 5ac313260ca39ce0b300bc779907f64f8fc49f0a Mon Sep 17 00:00:00 2001 From: guillaume Date: Fri, 14 Jul 2017 15:01:08 -0400 Subject: [PATCH 41/54] changed test to PythonAPItestharness and added tests on surface filter behavior at boundaries --- src/tracking.F90 | 1 - .../test_surface_tally/test_surface_tally.py | 156 +++++++++++++++++- 2 files changed, 154 insertions(+), 3 deletions(-) diff --git a/src/tracking.F90 b/src/tracking.F90 index 042e63c27f..546ae38809 100644 --- a/src/tracking.F90 +++ b/src/tracking.F90 @@ -36,7 +36,6 @@ contains integer :: next_level ! next coordinate level to check integer :: surface_crossed ! surface which particle is on integer :: lattice_translation(3) ! in-lattice translation vector - integer :: last_cell ! most recent cell particle was in integer :: n_event ! number of collisions/crossings real(8) :: d_boundary ! distance to nearest boundary real(8) :: d_collision ! sampled distance to collision diff --git a/tests/test_surface_tally/test_surface_tally.py b/tests/test_surface_tally/test_surface_tally.py index b04fcc6eba..d4d1c8aaa9 100644 --- a/tests/test_surface_tally/test_surface_tally.py +++ b/tests/test_surface_tally/test_surface_tally.py @@ -3,9 +3,161 @@ import os import sys sys.path.insert(0, os.pardir) -from testing_harness import TestHarness +from testing_harness import PyAPITestHarness +import openmc + + +class CreateSurfaceTallyTestHarness(PyAPITestHarness): + def _build_inputs(self): + # Instantiate some Materials and register the appropriate Nuclides + uo2 = openmc.Material(name='UO2 fuel at 2.4% wt enrichment') + uo2.set_density('g/cm3', 10.29769) + uo2.add_element('U', 1., enrichment=2.4) + uo2.add_element('O', 2.) + + borated_water = openmc.Material(name='Borated water') + borated_water.set_density('g/cm3', 0.740582) + borated_water.add_element('B', 4.0e-5) + borated_water.add_element('H', 5.0e-2) + borated_water.add_element('O', 2.4e-2) + borated_water.add_s_alpha_beta('c_H_in_H2O') + + # Instantiate a Materials collection and export to XML + materials_file = openmc.Materials([uo2, borated_water]) + + materials_file.export_to_xml() + + # Instantiate ZCylinder surfaces + fuel_or = openmc.ZCylinder(surface_id=1, x0=0, y0=0, R=0.4, name='Fuel OR') + left = openmc.XPlane(surface_id=2, x0=-0.62992, name='left') + right = openmc.XPlane(x0=0.62992, name='right') + bottom = openmc.YPlane(y0=-0.62992, name='bottom') + top = openmc.YPlane(y0=0.62992, name='top') + + left.boundary_type = 'vacuum' + right.boundary_type = 'reflective' + top.boundary_type = 'reflective' + bottom.boundary_type = 'reflective' + + # Instantiate Cells + fuel = openmc.Cell(name='fuel') + water = openmc.Cell(name='water') + + # Use surface half-spaces to define regions + fuel.region = -fuel_or + water.region = +fuel_or & -right & +bottom & -top + + # Register Materials with Cells + fuel.fill = uo2 + water.fill = borated_water + + # Instantiate pin cell Universe + pin_cell = openmc.Universe(name='pin cell') + pin_cell.add_cells([fuel, water]) + + # Instantiate root Cell and Universe + root_cell = openmc.Cell(name='root cell') + root_cell.region = +left & -right & +bottom & -top + root_cell.fill = pin_cell + root_univ = openmc.Universe(universe_id=0, name='root universe') + root_univ.add_cell(root_cell) + + # Instantiate a Geometry, register the root Universe, and export to XML + geometry = openmc.Geometry(root_univ) + geometry.export_to_xml() + + # Instantiate a Settings object, set all runtime parameters, and export to XML + settings_file = openmc.Settings() + settings_file.batches = 10 + settings_file.inactive = 5 + settings_file.particles = 1000 + + # Create an initial uniform spatial source distribution over fissionable zones + bounds = [-0.62992, -0.62992, -1, 0.62992, 0.62992, 1] + uniform_dist = openmc.stats.Box(bounds[:3], bounds[3:], only_fissionable=True) + settings_file.source = openmc.source.Source(space=uniform_dist) + settings_file.export_to_xml() + + # Tallies file + tallies_file = openmc.Tallies() + + # Cell to cell tallies + # These filters are same for all tallies + energy_filter = openmc.EnergyFilter(two_groups) + polar_filter = openmc.PolarFilter([0, np.pi / 4, np.pi]) + azimuthal_filter = openmc.AzimuthalFilter([0, np.pi / 4, np.pi]) + + + cell_to_cell_tallies = [] + tally_index = 0 + + for cell1 in pin_cell.get_all_cells(): + for cell2 in pin_cell.get_all_cells(): + + if cell1 != cell2 and abs(abs(cell1-cell2)-1) < 0.1: # no need for cell1 to cell1 + + # Cell to cell filters for partial current + cell_from_filter = openmc.CellFromFilter(cell1) + cell_to_filter = openmc.CellFilter(cell2) + + cell_to_cell_tallies.append(openmc.Tally(tally_id=2*tally_index, name=str(cell1)+'-'+str(cell2))) + cell_to_cell_tallies[2*tally_index].filters = [cell_from_filter, cell_to_filter, energy_filter, polar_filter, azimuthal_filter] + cell_to_cell_tallies[2*tally_index].scores = ['current'] + cell_to_cell_tallies[2*tally_index].estimator = 'analog' + + tallies_file.append(cell_to_cell_tallies[2*tally_index]) + + # Cell from + surface filters for partial current + surface_filter = openmc.SurfaceFilter([1]) + cell_to_cell_tallies.append(openmc.Tally(tally_id=2*tally_index+1, name=str(cell1)+'-surface1')) + cell_to_cell_tallies[2*tally_index+1].filters = [cell_from_filter, surface_filter, energy_filter, polar_filter, azimuthal_filter] + cell_to_cell_tallies[2*tally_index+1].scores = ['current'] + cell_to_cell_tallies[2*tally_index+1].estimator = 'analog' + tallies_file.append(cell_to_cell_tallies[2*tally_index+1]) + + tally_index += 1 + + # Surface filter on inner surface, for net current + surface_filter = openmc.SurfaceFilter([1]) + cell_to_cell_tallies.append(openmc.Tally(tally_id=2*tally_index, name='surface1')) + cell_to_cell_tallies[2*tally_index].filters = [surface_filter, energy_filter, polar_filter, azimuthal_filter] + cell_to_cell_tallies[2*tally_index].scores = ['current'] + cell_to_cell_tallies[2*tally_index].estimator = 'analog' + tallies_file.append(cell_to_cell_tallies[2*tally_index]) + + # Surface filter on left surface, vacuum BC, for net current = leakage + surface_filter = openmc.SurfaceFilter([1]) + cell_to_cell_tallies.append(openmc.Tally(tally_id=2*tally_index, name='surface1')) + cell_to_cell_tallies[2*tally_index].filters = [surface_filter, energy_filter, polar_filter, azimuthal_filter] + cell_to_cell_tallies[2*tally_index].scores = ['current'] + cell_to_cell_tallies[2*tally_index].estimator = 'analog' + tallies_file.append(cell_to_cell_tallies[2*tally_index]) + + # Surface filter on right surface, reflective, for net current = 0 + surface_filter = openmc.SurfaceFilter([1]) + cell_to_cell_tallies.append(openmc.Tally(tally_id=2*tally_index, name='surface1')) + cell_to_cell_tallies[2*tally_index].filters = [surface_filter, energy_filter, polar_filter, azimuthal_filter] + cell_to_cell_tallies[2*tally_index].scores = ['current'] + cell_to_cell_tallies[2*tally_index].estimator = 'analog' + tallies_file.append(cell_to_cell_tallies[2*tally_index]) + + tallies_file.export_to_xml() + + def _get_results(self): + """Digest info in the statepoint and return as a string.""" + # Read the statepoint file. + sp = openmc.StatePoint(self._sp_name) + + # Write out tally data. + outstr = '' + t = sp.get_tally() + outstr += 'tally {}:\n'.format(t.id) + outstr += 'sum = {:12.6E}\n'.format(t.sum[0, 0, 0]) + outstr += 'sum_sq = {:12.6E}\n'.format(t.sum_sq[0, 0, 0]) + + return outstr if __name__ == '__main__': - harness = TestHarness('statepoint.10.h5') + harness = CreateSurfaceTallyTestHarness('statepoint.10.h5') harness.main() From 87a1139f6232cc7521991c7f59e6cf69b5c85eb7 Mon Sep 17 00:00:00 2001 From: guillaume Date: Fri, 14 Jul 2017 15:13:27 -0400 Subject: [PATCH 42/54] cleaned up test and generated results --- tests/test_surface_tally/geometry.xml | 11 --- tests/test_surface_tally/materials.xml | 21 ----- tests/test_surface_tally/results_true.dat | 90 +------------------ tests/test_surface_tally/settings.xml | 12 --- tests/test_surface_tally/tallies.xml | 52 ----------- .../test_surface_tally/test_surface_tally.py | 59 ++++++------ 6 files changed, 32 insertions(+), 213 deletions(-) delete mode 100644 tests/test_surface_tally/geometry.xml delete mode 100644 tests/test_surface_tally/materials.xml delete mode 100644 tests/test_surface_tally/settings.xml delete mode 100644 tests/test_surface_tally/tallies.xml diff --git a/tests/test_surface_tally/geometry.xml b/tests/test_surface_tally/geometry.xml deleted file mode 100644 index 19eae6333b..0000000000 --- a/tests/test_surface_tally/geometry.xml +++ /dev/null @@ -1,11 +0,0 @@ - - - - - - - - - - - diff --git a/tests/test_surface_tally/materials.xml b/tests/test_surface_tally/materials.xml deleted file mode 100644 index 506f5f74ae..0000000000 --- a/tests/test_surface_tally/materials.xml +++ /dev/null @@ -1,21 +0,0 @@ - - - - - - - - - - - - - - - - - - - - - diff --git a/tests/test_surface_tally/results_true.dat b/tests/test_surface_tally/results_true.dat index f9f473b19c..e47a0635c8 100644 --- a/tests/test_surface_tally/results_true.dat +++ b/tests/test_surface_tally/results_true.dat @@ -1,87 +1,3 @@ -k-combined: -1.177844E+00 1.713522E-02 -tally 1: -9.400000E-01 -1.778300E-01 -2.561000E+00 -1.315759E+00 -8.618000E+00 -1.487929E+01 -2.559600E+01 -1.311144E+02 -4.700000E-02 -4.870000E-04 -1.010000E-01 -2.207000E-03 -3.780000E-01 -2.889000E-02 -1.089000E+00 -2.432570E-01 -tally 2: -9.400000E-01 -1.778300E-01 -2.561000E+00 -1.315759E+00 -8.618000E+00 -1.487929E+01 -2.559600E+01 -1.311144E+02 -4.700000E-02 -4.870000E-04 -1.010000E-01 -2.207000E-03 -3.780000E-01 -2.889000E-02 -1.089000E+00 -2.432570E-01 -tally 3: -9.070000E-01 -1.648450E-01 -2.495000E+00 -1.247251E+00 -8.515000E+00 -1.452994E+01 -2.539500E+01 -1.290580E+02 -3.800000E-02 -3.220000E-04 -8.400000E-02 -1.480000E-03 -3.540000E-01 -2.528400E-02 -9.770000E-01 -1.957910E-01 -tally 4: -9.070000E-01 -1.648450E-01 -2.495000E+00 -1.247251E+00 -8.515000E+00 -1.452994E+01 -2.539500E+01 -1.290580E+02 -3.800000E-02 -3.220000E-04 -8.400000E-02 -1.480000E-03 -3.540000E-01 -2.528400E-02 -9.770000E-01 -1.957910E-01 -tally 5: -1.847000E+00 -6.843910E-01 -5.056000E+00 -5.124718E+00 -1.713300E+01 -5.881492E+01 -5.099100E+01 -5.203348E+02 -8.500000E-02 -1.595000E-03 -1.850000E-01 -7.259000E-03 -7.320000E-01 -1.080440E-01 -2.066000E+00 -8.753660E-01 +tally 0: +sum = 8.800000E-02 +sum_sq = 1.582000E-03 diff --git a/tests/test_surface_tally/settings.xml b/tests/test_surface_tally/settings.xml deleted file mode 100644 index d97ee454a4..0000000000 --- a/tests/test_surface_tally/settings.xml +++ /dev/null @@ -1,12 +0,0 @@ - - - eigenvalue - 1000 - 10 - 5 - - - -0.62992 -0.62992 -1 0.62992 0.62992 1 - - - diff --git a/tests/test_surface_tally/tallies.xml b/tests/test_surface_tally/tallies.xml deleted file mode 100644 index fd76a84d22..0000000000 --- a/tests/test_surface_tally/tallies.xml +++ /dev/null @@ -1,52 +0,0 @@ - - - - 10000 - - - 10001 - - - 0.0 4000000.0 20000000.0 - - - 0.0 0.785398163397 3.14159265359 - - - 0.0 0.785398163397 3.14159265359 - - - 1 - - - 10001 - - - 10000 - - - 10003 10004 10000 10001 10002 - total-current - analog - - - 10003 10005 10000 10001 10002 - total-current - analog - - - 10006 10007 10000 10001 10002 - total-current - analog - - - 10006 10005 10000 10001 10002 - total-current - analog - - - 10005 10000 10001 10002 - total-current - analog - - diff --git a/tests/test_surface_tally/test_surface_tally.py b/tests/test_surface_tally/test_surface_tally.py index d4d1c8aaa9..31850462ce 100644 --- a/tests/test_surface_tally/test_surface_tally.py +++ b/tests/test_surface_tally/test_surface_tally.py @@ -4,6 +4,7 @@ import os import sys sys.path.insert(0, os.pardir) from testing_harness import PyAPITestHarness +import numpy as np import openmc @@ -30,7 +31,7 @@ class CreateSurfaceTallyTestHarness(PyAPITestHarness): # Instantiate ZCylinder surfaces fuel_or = openmc.ZCylinder(surface_id=1, x0=0, y0=0, R=0.4, name='Fuel OR') left = openmc.XPlane(surface_id=2, x0=-0.62992, name='left') - right = openmc.XPlane(x0=0.62992, name='right') + right = openmc.XPlane(surface_id=3, x0=0.62992, name='right') bottom = openmc.YPlane(y0=-0.62992, name='bottom') top = openmc.YPlane(y0=0.62992, name='top') @@ -83,6 +84,7 @@ class CreateSurfaceTallyTestHarness(PyAPITestHarness): # Cell to cell tallies # These filters are same for all tallies + two_groups = np.array([0., 4, 20.]) * 1e6 energy_filter = openmc.EnergyFilter(two_groups) polar_filter = openmc.PolarFilter([0, np.pi / 4, np.pi]) azimuthal_filter = openmc.AzimuthalFilter([0, np.pi / 4, np.pi]) @@ -100,46 +102,43 @@ class CreateSurfaceTallyTestHarness(PyAPITestHarness): cell_from_filter = openmc.CellFromFilter(cell1) cell_to_filter = openmc.CellFilter(cell2) - cell_to_cell_tallies.append(openmc.Tally(tally_id=2*tally_index, name=str(cell1)+'-'+str(cell2))) - cell_to_cell_tallies[2*tally_index].filters = [cell_from_filter, cell_to_filter, energy_filter, polar_filter, azimuthal_filter] - cell_to_cell_tallies[2*tally_index].scores = ['current'] - cell_to_cell_tallies[2*tally_index].estimator = 'analog' - - tallies_file.append(cell_to_cell_tallies[2*tally_index]) + cell_to_cell_tallies.append(openmc.Tally(tally_id=tally_index, name=str(cell1)+'-'+str(cell2))) + cell_to_cell_tallies[tally_index].filters = [cell_from_filter, cell_to_filter, energy_filter, polar_filter, azimuthal_filter] + cell_to_cell_tallies[tally_index].scores = ['current'] + tallies_file.append(cell_to_cell_tallies[tally_index]) + tally_index += 1 # Cell from + surface filters for partial current surface_filter = openmc.SurfaceFilter([1]) - cell_to_cell_tallies.append(openmc.Tally(tally_id=2*tally_index+1, name=str(cell1)+'-surface1')) - cell_to_cell_tallies[2*tally_index+1].filters = [cell_from_filter, surface_filter, energy_filter, polar_filter, azimuthal_filter] - cell_to_cell_tallies[2*tally_index+1].scores = ['current'] - cell_to_cell_tallies[2*tally_index+1].estimator = 'analog' - tallies_file.append(cell_to_cell_tallies[2*tally_index+1]) - + cell_to_cell_tallies.append(openmc.Tally(tally_id=tally_index, name=str(cell1)+'-surface1')) + cell_to_cell_tallies[tally_index].filters = [cell_from_filter, surface_filter, energy_filter, polar_filter, azimuthal_filter] + cell_to_cell_tallies[tally_index].scores = ['current'] + tallies_file.append(cell_to_cell_tallies[tally_index]) tally_index += 1 # Surface filter on inner surface, for net current surface_filter = openmc.SurfaceFilter([1]) - cell_to_cell_tallies.append(openmc.Tally(tally_id=2*tally_index, name='surface1')) - cell_to_cell_tallies[2*tally_index].filters = [surface_filter, energy_filter, polar_filter, azimuthal_filter] - cell_to_cell_tallies[2*tally_index].scores = ['current'] - cell_to_cell_tallies[2*tally_index].estimator = 'analog' - tallies_file.append(cell_to_cell_tallies[2*tally_index]) + surf_tally1 = openmc.Tally(tally_id=tally_index, name='net_cylinder') + surf_tally1.filters = [surface_filter, energy_filter, polar_filter, azimuthal_filter] + surf_tally1.scores = ['current'] + tallies_file.append(surf_tally1) + tally_index += 1 # Surface filter on left surface, vacuum BC, for net current = leakage - surface_filter = openmc.SurfaceFilter([1]) - cell_to_cell_tallies.append(openmc.Tally(tally_id=2*tally_index, name='surface1')) - cell_to_cell_tallies[2*tally_index].filters = [surface_filter, energy_filter, polar_filter, azimuthal_filter] - cell_to_cell_tallies[2*tally_index].scores = ['current'] - cell_to_cell_tallies[2*tally_index].estimator = 'analog' - tallies_file.append(cell_to_cell_tallies[2*tally_index]) + surface_filter = openmc.SurfaceFilter([2]) + surf_tally2 = openmc.Tally(tally_id=tally_index, name='leakage_left') + surf_tally2.filters = [surface_filter, energy_filter, polar_filter, azimuthal_filter] + surf_tally2.scores = ['current'] + tallies_file.append(surf_tally2) + tally_index += 1 # Surface filter on right surface, reflective, for net current = 0 - surface_filter = openmc.SurfaceFilter([1]) - cell_to_cell_tallies.append(openmc.Tally(tally_id=2*tally_index, name='surface1')) - cell_to_cell_tallies[2*tally_index].filters = [surface_filter, energy_filter, polar_filter, azimuthal_filter] - cell_to_cell_tallies[2*tally_index].scores = ['current'] - cell_to_cell_tallies[2*tally_index].estimator = 'analog' - tallies_file.append(cell_to_cell_tallies[2*tally_index]) + surface_filter = openmc.SurfaceFilter([3]) + surf_tally3 = openmc.Tally(tally_id=tally_index, name='net_right') + surf_tally3.filters = [surface_filter, energy_filter, polar_filter, azimuthal_filter] + surf_tally3.scores = ['current'] + tallies_file.append(surf_tally3) + tally_index += 1 tallies_file.export_to_xml() From 7a2ce111c6a580d3f1b2ba46732e611a4d885a91 Mon Sep 17 00:00:00 2001 From: guillaume Date: Fri, 14 Jul 2017 15:14:38 -0400 Subject: [PATCH 43/54] added inputs_true to git --- tests/test_surface_tally/inputs_true.dat | 105 +++++++++++++++++++++++ 1 file changed, 105 insertions(+) create mode 100644 tests/test_surface_tally/inputs_true.dat diff --git a/tests/test_surface_tally/inputs_true.dat b/tests/test_surface_tally/inputs_true.dat new file mode 100644 index 0000000000..b1a52992d3 --- /dev/null +++ b/tests/test_surface_tally/inputs_true.dat @@ -0,0 +1,105 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + eigenvalue + 1000 + 10 + 5 + + + -0.62992 -0.62992 -1 0.62992 0.62992 1 + + + + + + + 10000 + + + 10001 + + + 0.0 4000000.0 20000000.0 + + + 0.0 0.785398163397 3.14159265359 + + + 0.0 0.785398163397 3.14159265359 + + + 1 + + + 10001 + + + 10000 + + + 2 + + + 3 + + + 10003 10004 10000 10001 10002 + current + + + 10003 10005 10000 10001 10002 + current + + + 10006 10007 10000 10001 10002 + current + + + 10006 10005 10000 10001 10002 + current + + + 10005 10000 10001 10002 + current + + + 10010 10000 10001 10002 + current + + + 10011 10000 10001 10002 + current + + From 74f56063582ee380b23ad97c872283f8c4efbaeb Mon Sep 17 00:00:00 2001 From: guillaume Date: Fri, 14 Jul 2017 15:30:02 -0400 Subject: [PATCH 44/54] white space 27 lines before what check_source was saying --- src/input_xml.F90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/input_xml.F90 b/src/input_xml.F90 index e0fa63e0db..dfc2e853c5 100644 --- a/src/input_xml.F90 +++ b/src/input_xml.F90 @@ -3810,7 +3810,7 @@ contains t % score_bins(j) = SCORE_FISS_Q_RECOV case ('current') - ! Check which type of current is desired: mesh currents or + ! Check which type of current is desired: mesh currents or ! surface currents if (t % find_filter(FILTER_SURFACE) > 0 .or. & &t % find_filter(FILTER_CELL) > 0 .or. & From d28bb4f9b38a99f0686b21b762428bc151026aac Mon Sep 17 00:00:00 2001 From: guillaume Date: Fri, 14 Jul 2017 15:59:28 -0400 Subject: [PATCH 45/54] get_results now does all tallies instead of 1 --- tests/test_surface_tally/inputs_true.dat | 4 +- tests/test_surface_tally/results_true.dat | 60 ++++++++++++++++++- .../test_surface_tally/test_surface_tally.py | 20 ++++--- 3 files changed, 70 insertions(+), 14 deletions(-) diff --git a/tests/test_surface_tally/inputs_true.dat b/tests/test_surface_tally/inputs_true.dat index b1a52992d3..d161c88dc1 100644 --- a/tests/test_surface_tally/inputs_true.dat +++ b/tests/test_surface_tally/inputs_true.dat @@ -33,9 +33,9 @@ eigenvalue - 1000 + 100000 10 - 5 + 0 -0.62992 -0.62992 -1 0.62992 0.62992 1 diff --git a/tests/test_surface_tally/results_true.dat b/tests/test_surface_tally/results_true.dat index e47a0635c8..cb9e764757 100644 --- a/tests/test_surface_tally/results_true.dat +++ b/tests/test_surface_tally/results_true.dat @@ -1,3 +1,57 @@ -tally 0: -sum = 8.800000E-02 -sum_sq = 1.582000E-03 +d_material,d_nuclide,d_variable,score,mean,std. dev. +,,,current,2.0521000e-02,1.8594772e-04 +,,,current,9.4192000e-02,4.2316217e-04 +,,,current,1.2748000e-01,5.1128596e-04 +,,,current,7.2691500e-01,2.2061909e-03 +,,,current,1.8240000e-03,4.3389707e-05 +,,,current,9.6750000e-03,1.7794662e-04 +,,,current,1.1288000e-02,2.4212394e-04 +,,,current,6.9557000e-02,8.3536692e-04 +,,,current,2.0521000e-02,1.8594772e-04 +,,,current,9.4192000e-02,4.2316217e-04 +,,,current,1.2748000e-01,5.1128596e-04 +,,,current,7.2691500e-01,2.2061909e-03 +,,,current,1.8240000e-03,4.3389707e-05 +,,,current,9.6750000e-03,1.7794662e-04 +,,,current,1.1288000e-02,2.4212394e-04 +,,,current,6.9557000e-02,8.3536692e-04 +,,,current,2.4030000e-03,6.0443730e-05 +,,,current,5.4706000e-02,2.8226150e-04 +,,,current,1.0554000e-02,9.9712921e-05 +,,,current,4.5515300e-01,1.1105905e-03 +,,,current,2.1000000e-05,5.0442487e-06 +,,,current,5.2200000e-03,1.2808851e-04 +,,,current,9.1000000e-05,9.0000000e-06 +,,,current,3.9997000e-02,4.9880981e-04 +,,,current,-2.4030000e-03,6.0443730e-05 +,,,current,-5.4706000e-02,2.8226150e-04 +,,,current,-1.0554000e-02,9.9712921e-05 +,,,current,-4.5515300e-01,1.1105905e-03 +,,,current,-2.1000000e-05,5.0442487e-06 +,,,current,-5.2200000e-03,1.2808851e-04 +,,,current,-9.1000000e-05,9.0000000e-06 +,,,current,-3.9997000e-02,4.9880981e-04 +,,,current,1.8118000e-02,1.6043552e-04 +,,,current,3.9486000e-02,2.7664337e-04 +,,,current,1.1692600e-01,4.8142889e-04 +,,,current,2.7176200e-01,1.5048527e-03 +,,,current,1.8030000e-03,4.0525712e-05 +,,,current,4.4550000e-03,1.6878816e-04 +,,,current,1.1197000e-02,2.4260187e-04 +,,,current,2.9560000e-02,4.5309307e-04 +,,,current,0.0000000e+00,0.0000000e+00 +,,,current,-4.9065000e-02,3.1853222e-04 +,,,current,0.0000000e+00,0.0000000e+00 +,,,current,-4.0249600e-01,3.6746338e-04 +,,,current,0.0000000e+00,0.0000000e+00 +,,,current,-4.9020000e-03,7.5627743e-05 +,,,current,0.0000000e+00,0.0000000e+00 +,,,current,-3.5016000e-02,2.8455501e-04 +,,,current,0.0000000e+00,0.0000000e+00 +,,,current,-3.5893000e-02,1.9351170e-04 +,,,current,0.0000000e+00,0.0000000e+00 +,,,current,-2.5851100e-01,7.8233617e-04 +,,,current,0.0000000e+00,0.0000000e+00 +,,,current,-2.9840000e-03,7.1277237e-05 +,,,current,0.0000000e+00,0.0000000e+00 +,,,current,-2.0319000e-02,2.5016195e-04 diff --git a/tests/test_surface_tally/test_surface_tally.py b/tests/test_surface_tally/test_surface_tally.py index 31850462ce..dfed786151 100644 --- a/tests/test_surface_tally/test_surface_tally.py +++ b/tests/test_surface_tally/test_surface_tally.py @@ -6,7 +6,7 @@ sys.path.insert(0, os.pardir) from testing_harness import PyAPITestHarness import numpy as np import openmc - +import pandas as pd class CreateSurfaceTallyTestHarness(PyAPITestHarness): def _build_inputs(self): @@ -70,8 +70,8 @@ class CreateSurfaceTallyTestHarness(PyAPITestHarness): # Instantiate a Settings object, set all runtime parameters, and export to XML settings_file = openmc.Settings() settings_file.batches = 10 - settings_file.inactive = 5 - settings_file.particles = 1000 + settings_file.inactive = 0 + settings_file.particles = 100000 # Create an initial uniform spatial source distribution over fissionable zones bounds = [-0.62992, -0.62992, -1, 0.62992, 0.62992, 1] @@ -147,13 +147,15 @@ class CreateSurfaceTallyTestHarness(PyAPITestHarness): # Read the statepoint file. sp = openmc.StatePoint(self._sp_name) - # Write out tally data. - outstr = '' - t = sp.get_tally() - outstr += 'tally {}:\n'.format(t.id) - outstr += 'sum = {:12.6E}\n'.format(t.sum[0, 0, 0]) - outstr += 'sum_sq = {:12.6E}\n'.format(t.sum_sq[0, 0, 0]) + # Extract the tally data as a Pandas DataFrame. + df = pd.DataFrame() + for t in sp.tallies.values(): + df = df.append(t.get_pandas_dataframe(), ignore_index=True) + # Extract the relevant data as a CSV string. + cols = ('d_material', 'd_nuclide', 'd_variable', 'score', 'mean', + 'std. dev.') + return df.to_csv(None, columns=cols, index=False, float_format='%.7e') return outstr From 1b4f2ae0254c1f364b0d15a2c80585772dd83545 Mon Sep 17 00:00:00 2001 From: guillaume Date: Fri, 14 Jul 2017 16:19:58 -0400 Subject: [PATCH 46/54] not much --- tests/test_surface_tally/inputs_true.dat | 2 +- tests/test_surface_tally/results_true.dat | 96 +++++++++---------- .../test_surface_tally/test_surface_tally.py | 2 +- 3 files changed, 50 insertions(+), 50 deletions(-) diff --git a/tests/test_surface_tally/inputs_true.dat b/tests/test_surface_tally/inputs_true.dat index d161c88dc1..28c216075e 100644 --- a/tests/test_surface_tally/inputs_true.dat +++ b/tests/test_surface_tally/inputs_true.dat @@ -33,7 +33,7 @@ eigenvalue - 100000 + 1000 10 0 diff --git a/tests/test_surface_tally/results_true.dat b/tests/test_surface_tally/results_true.dat index cb9e764757..41e5484756 100644 --- a/tests/test_surface_tally/results_true.dat +++ b/tests/test_surface_tally/results_true.dat @@ -1,57 +1,57 @@ d_material,d_nuclide,d_variable,score,mean,std. dev. -,,,current,2.0521000e-02,1.8594772e-04 -,,,current,9.4192000e-02,4.2316217e-04 -,,,current,1.2748000e-01,5.1128596e-04 -,,,current,7.2691500e-01,2.2061909e-03 -,,,current,1.8240000e-03,4.3389707e-05 -,,,current,9.6750000e-03,1.7794662e-04 -,,,current,1.1288000e-02,2.4212394e-04 -,,,current,6.9557000e-02,8.3536692e-04 -,,,current,2.0521000e-02,1.8594772e-04 -,,,current,9.4192000e-02,4.2316217e-04 -,,,current,1.2748000e-01,5.1128596e-04 -,,,current,7.2691500e-01,2.2061909e-03 -,,,current,1.8240000e-03,4.3389707e-05 -,,,current,9.6750000e-03,1.7794662e-04 -,,,current,1.1288000e-02,2.4212394e-04 -,,,current,6.9557000e-02,8.3536692e-04 -,,,current,2.4030000e-03,6.0443730e-05 -,,,current,5.4706000e-02,2.8226150e-04 -,,,current,1.0554000e-02,9.9712921e-05 -,,,current,4.5515300e-01,1.1105905e-03 -,,,current,2.1000000e-05,5.0442487e-06 -,,,current,5.2200000e-03,1.2808851e-04 -,,,current,9.1000000e-05,9.0000000e-06 -,,,current,3.9997000e-02,4.9880981e-04 -,,,current,-2.4030000e-03,6.0443730e-05 -,,,current,-5.4706000e-02,2.8226150e-04 -,,,current,-1.0554000e-02,9.9712921e-05 -,,,current,-4.5515300e-01,1.1105905e-03 -,,,current,-2.1000000e-05,5.0442487e-06 -,,,current,-5.2200000e-03,1.2808851e-04 -,,,current,-9.1000000e-05,9.0000000e-06 -,,,current,-3.9997000e-02,4.9880981e-04 -,,,current,1.8118000e-02,1.6043552e-04 -,,,current,3.9486000e-02,2.7664337e-04 -,,,current,1.1692600e-01,4.8142889e-04 -,,,current,2.7176200e-01,1.5048527e-03 -,,,current,1.8030000e-03,4.0525712e-05 -,,,current,4.4550000e-03,1.6878816e-04 -,,,current,1.1197000e-02,2.4260187e-04 -,,,current,2.9560000e-02,4.5309307e-04 +,,,current,2.4000000e-02,5.6725460e-03 +,,,current,1.0080000e-01,7.0266161e-03 +,,,current,1.2900000e-01,3.1304952e-03 +,,,current,7.2540000e-01,1.4240943e-02 +,,,current,1.8000000e-03,4.1633320e-04 +,,,current,8.3000000e-03,1.4146063e-03 +,,,current,1.5300000e-02,2.9666667e-03 +,,,current,7.9700000e-02,4.3717528e-03 +,,,current,2.4000000e-02,5.6725460e-03 +,,,current,1.0080000e-01,7.0266161e-03 +,,,current,1.2900000e-01,3.1304952e-03 +,,,current,7.2540000e-01,1.4240943e-02 +,,,current,1.8000000e-03,4.1633320e-04 +,,,current,8.3000000e-03,1.4146063e-03 +,,,current,1.5300000e-02,2.9666667e-03 +,,,current,7.9700000e-02,4.3717528e-03 +,,,current,1.5000000e-03,4.0138649e-04 +,,,current,5.7500000e-02,4.1075811e-03 +,,,current,1.0000000e-02,1.3416408e-03 +,,,current,4.5890000e-01,8.6941232e-03 ,,,current,0.0000000e+00,0.0000000e+00 -,,,current,-4.9065000e-02,3.1853222e-04 +,,,current,4.6000000e-03,1.0241528e-03 +,,,current,2.0000000e-04,1.3333333e-04 +,,,current,4.6700000e-02,3.5654203e-03 +,,,current,-1.5000000e-03,4.0138649e-04 +,,,current,-5.7500000e-02,4.1075811e-03 +,,,current,-1.0000000e-02,1.3416408e-03 +,,,current,-4.5890000e-01,8.6941232e-03 ,,,current,0.0000000e+00,0.0000000e+00 -,,,current,-4.0249600e-01,3.6746338e-04 +,,,current,-4.6000000e-03,1.0241528e-03 +,,,current,-2.0000000e-04,1.3333333e-04 +,,,current,-4.6700000e-02,3.5654203e-03 +,,,current,2.2500000e-02,5.6651959e-03 +,,,current,4.3300000e-02,4.1714639e-03 +,,,current,1.1900000e-01,2.8635642e-03 +,,,current,2.6650000e-01,1.4426249e-02 +,,,current,1.8000000e-03,4.1633320e-04 +,,,current,3.7000000e-03,1.0333333e-03 +,,,current,1.5100000e-02,2.8691849e-03 +,,,current,3.3000000e-02,1.7448336e-03 ,,,current,0.0000000e+00,0.0000000e+00 -,,,current,-4.9020000e-03,7.5627743e-05 +,,,current,-5.2600000e-02,4.8355167e-03 ,,,current,0.0000000e+00,0.0000000e+00 -,,,current,-3.5016000e-02,2.8455501e-04 +,,,current,-3.9890000e-01,8.1150751e-03 ,,,current,0.0000000e+00,0.0000000e+00 -,,,current,-3.5893000e-02,1.9351170e-04 +,,,current,-4.2000000e-03,7.1180522e-04 ,,,current,0.0000000e+00,0.0000000e+00 -,,,current,-2.5851100e-01,7.8233617e-04 +,,,current,-3.9000000e-02,4.1686662e-03 ,,,current,0.0000000e+00,0.0000000e+00 -,,,current,-2.9840000e-03,7.1277237e-05 +,,,current,-4.0200000e-02,4.2656249e-03 ,,,current,0.0000000e+00,0.0000000e+00 -,,,current,-2.0319000e-02,2.5016195e-04 +,,,current,-2.6310000e-01,1.0013824e-02 +,,,current,0.0000000e+00,0.0000000e+00 +,,,current,-2.6000000e-03,4.0000000e-04 +,,,current,0.0000000e+00,0.0000000e+00 +,,,current,-2.6700000e-02,5.2218558e-03 diff --git a/tests/test_surface_tally/test_surface_tally.py b/tests/test_surface_tally/test_surface_tally.py index dfed786151..c8db320f93 100644 --- a/tests/test_surface_tally/test_surface_tally.py +++ b/tests/test_surface_tally/test_surface_tally.py @@ -71,7 +71,7 @@ class CreateSurfaceTallyTestHarness(PyAPITestHarness): settings_file = openmc.Settings() settings_file.batches = 10 settings_file.inactive = 0 - settings_file.particles = 100000 + settings_file.particles = 1000 # Create an initial uniform spatial source distribution over fissionable zones bounds = [-0.62992, -0.62992, -1, 0.62992, 0.62992, 1] From 937c5f6f5d68ebf51062929967f123b4800d7603 Mon Sep 17 00:00:00 2001 From: guillaume Date: Fri, 14 Jul 2017 16:50:59 -0400 Subject: [PATCH 47/54] rewording --- tests/test_surface_tally/test_surface_tally.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/tests/test_surface_tally/test_surface_tally.py b/tests/test_surface_tally/test_surface_tally.py index c8db320f93..caa0737c5b 100644 --- a/tests/test_surface_tally/test_surface_tally.py +++ b/tests/test_surface_tally/test_surface_tally.py @@ -132,10 +132,10 @@ class CreateSurfaceTallyTestHarness(PyAPITestHarness): tallies_file.append(surf_tally2) tally_index += 1 - # Surface filter on right surface, reflective, for net current = 0 + # Surface filter on right surface, reflective, current tally doesnt pick up the zero net current though surface_filter = openmc.SurfaceFilter([3]) surf_tally3 = openmc.Tally(tally_id=tally_index, name='net_right') - surf_tally3.filters = [surface_filter, energy_filter, polar_filter, azimuthal_filter] + surf_tally3.filters = [surface_filter, energy_filterr] surf_tally3.scores = ['current'] tallies_file.append(surf_tally3) tally_index += 1 From 5b2aae3b05b29fa9650a51fd54d47395d7e56b6c Mon Sep 17 00:00:00 2001 From: guillaume Date: Fri, 14 Jul 2017 17:17:42 -0400 Subject: [PATCH 48/54] debugging test --- tests/test_surface_tally/inputs_true.dat | 2 +- tests/test_surface_tally/results_true.dat | 10 ++-------- tests/test_surface_tally/test_surface_tally.py | 2 +- 3 files changed, 4 insertions(+), 10 deletions(-) diff --git a/tests/test_surface_tally/inputs_true.dat b/tests/test_surface_tally/inputs_true.dat index 28c216075e..213097e8d3 100644 --- a/tests/test_surface_tally/inputs_true.dat +++ b/tests/test_surface_tally/inputs_true.dat @@ -99,7 +99,7 @@ current - 10011 10000 10001 10002 + 10011 10000 current diff --git a/tests/test_surface_tally/results_true.dat b/tests/test_surface_tally/results_true.dat index 41e5484756..df8c1597f4 100644 --- a/tests/test_surface_tally/results_true.dat +++ b/tests/test_surface_tally/results_true.dat @@ -47,11 +47,5 @@ d_material,d_nuclide,d_variable,score,mean,std. dev. ,,,current,-4.2000000e-03,7.1180522e-04 ,,,current,0.0000000e+00,0.0000000e+00 ,,,current,-3.9000000e-02,4.1686662e-03 -,,,current,0.0000000e+00,0.0000000e+00 -,,,current,-4.0200000e-02,4.2656249e-03 -,,,current,0.0000000e+00,0.0000000e+00 -,,,current,-2.6310000e-01,1.0013824e-02 -,,,current,0.0000000e+00,0.0000000e+00 -,,,current,-2.6000000e-03,4.0000000e-04 -,,,current,0.0000000e+00,0.0000000e+00 -,,,current,-2.6700000e-02,5.2218558e-03 +,,,current,-5.9230000e-01,1.2109730e-02 +,,,current,-5.1900000e-02,3.3381632e-03 diff --git a/tests/test_surface_tally/test_surface_tally.py b/tests/test_surface_tally/test_surface_tally.py index caa0737c5b..cafdd84e23 100644 --- a/tests/test_surface_tally/test_surface_tally.py +++ b/tests/test_surface_tally/test_surface_tally.py @@ -135,7 +135,7 @@ class CreateSurfaceTallyTestHarness(PyAPITestHarness): # Surface filter on right surface, reflective, current tally doesnt pick up the zero net current though surface_filter = openmc.SurfaceFilter([3]) surf_tally3 = openmc.Tally(tally_id=tally_index, name='net_right') - surf_tally3.filters = [surface_filter, energy_filterr] + surf_tally3.filters = [surface_filter, energy_filter] surf_tally3.scores = ['current'] tallies_file.append(surf_tally3) tally_index += 1 From f0c051db5f5779b3b4a7cc7ba987c4eacbdfb276 Mon Sep 17 00:00:00 2001 From: guillaume Date: Sun, 6 Aug 2017 13:54:03 -0400 Subject: [PATCH 49/54] corrected syntax in PyAPI test, and deleted forgotten cellto options in input handling --- openmc/filter.py | 3 +- src/input_xml.F90 | 2 +- .../test_surface_tally/test_surface_tally.py | 145 ++++++++++-------- 3 files changed, 84 insertions(+), 66 deletions(-) diff --git a/openmc/filter.py b/openmc/filter.py index e822b9553e..514e96f261 100644 --- a/openmc/filter.py +++ b/openmc/filter.py @@ -18,8 +18,7 @@ AUTO_FILTER_ID = 10000 _FILTER_TYPES = ['universe', 'material', 'cell', 'cellborn', 'surface', 'mesh', 'energy', 'energyout', 'mu', 'polar', 'azimuthal', - 'distribcell', 'delayedgroup', 'energyfunction', 'cellfrom', - 'cellto'] + 'distribcell', 'delayedgroup', 'energyfunction', 'cellfrom'] _CURRENT_NAMES = {1: 'x-min out', 2: 'x-min in', 3: 'x-max out', 4: 'x-max in', diff --git a/src/input_xml.F90 b/src/input_xml.F90 index dfc2e853c5..cfb80b4fa5 100644 --- a/src/input_xml.F90 +++ b/src/input_xml.F90 @@ -3004,7 +3004,7 @@ contains end if n_words = node_word_count(node_filt, "bins") case ("mesh", "universe", "material", "cell", "distribcell", & - "cellborn", "cellto", "cellfrom", "surface", "delayedgroup") + "cellborn", "cellfrom", "surface", "delayedgroup") if (.not. check_for_node(node_filt, "bins")) then call fatal_error("Bins not set in filter " // trim(to_str(filter_id))) end if diff --git a/tests/test_surface_tally/test_surface_tally.py b/tests/test_surface_tally/test_surface_tally.py index cafdd84e23..5ebe30d587 100644 --- a/tests/test_surface_tally/test_surface_tally.py +++ b/tests/test_surface_tally/test_surface_tally.py @@ -8,7 +8,8 @@ import numpy as np import openmc import pandas as pd -class CreateSurfaceTallyTestHarness(PyAPITestHarness): + +class SurfaceTallyTestHarness(PyAPITestHarness): def _build_inputs(self): # Instantiate some Materials and register the appropriate Nuclides uo2 = openmc.Material(name='UO2 fuel at 2.4% wt enrichment') @@ -25,15 +26,15 @@ class CreateSurfaceTallyTestHarness(PyAPITestHarness): # Instantiate a Materials collection and export to XML materials_file = openmc.Materials([uo2, borated_water]) - materials_file.export_to_xml() # Instantiate ZCylinder surfaces - fuel_or = openmc.ZCylinder(surface_id=1, x0=0, y0=0, R=0.4, name='Fuel OR') - left = openmc.XPlane(surface_id=2, x0=-0.62992, name='left') - right = openmc.XPlane(surface_id=3, x0=0.62992, name='right') - bottom = openmc.YPlane(y0=-0.62992, name='bottom') - top = openmc.YPlane(y0=0.62992, name='top') + fuel_or = openmc.ZCylinder(surface_id=1, x0=0, y0=0, R=1, \ + name='Fuel OR') + left = openmc.XPlane(surface_id=2, x0=-2, name='left') + right = openmc.XPlane(surface_id=3, x0=2, name='right') + bottom = openmc.YPlane(y0=-2, name='bottom') + top = openmc.YPlane(y0=2, name='top') left.boundary_type = 'vacuum' right.boundary_type = 'reflective' @@ -49,7 +50,7 @@ class CreateSurfaceTallyTestHarness(PyAPITestHarness): water.region = +fuel_or & -right & +bottom & -top # Register Materials with Cells - fuel.fill = uo2 + fuel.fill = uo2 water.fill = borated_water # Instantiate pin cell Universe @@ -63,82 +64,101 @@ class CreateSurfaceTallyTestHarness(PyAPITestHarness): root_univ = openmc.Universe(universe_id=0, name='root universe') root_univ.add_cell(root_cell) - # Instantiate a Geometry, register the root Universe, and export to XML + # Instantiate a Geometry, register the root Universe geometry = openmc.Geometry(root_univ) geometry.export_to_xml() - # Instantiate a Settings object, set all runtime parameters, and export to XML + # Instantiate a Settings object, set all runtime parameters settings_file = openmc.Settings() settings_file.batches = 10 settings_file.inactive = 0 settings_file.particles = 1000 + #settings_file.output = {'tallies': True} - # Create an initial uniform spatial source distribution over fissionable zones + # Create an initial uniform spatial source distribution bounds = [-0.62992, -0.62992, -1, 0.62992, 0.62992, 1] - uniform_dist = openmc.stats.Box(bounds[:3], bounds[3:], only_fissionable=True) + uniform_dist = openmc.stats.Box(bounds[:3], bounds[3:],\ + only_fissionable=True) settings_file.source = openmc.source.Source(space=uniform_dist) settings_file.export_to_xml() # Tallies file tallies_file = openmc.Tallies() - # Cell to cell tallies - # These filters are same for all tallies + # Create partial current tallies from fuel to water + # Filters two_groups = np.array([0., 4, 20.]) * 1e6 - energy_filter = openmc.EnergyFilter(two_groups) - polar_filter = openmc.PolarFilter([0, np.pi / 4, np.pi]) - azimuthal_filter = openmc.AzimuthalFilter([0, np.pi / 4, np.pi]) - - - cell_to_cell_tallies = [] - tally_index = 0 - - for cell1 in pin_cell.get_all_cells(): - for cell2 in pin_cell.get_all_cells(): - - if cell1 != cell2 and abs(abs(cell1-cell2)-1) < 0.1: # no need for cell1 to cell1 - - # Cell to cell filters for partial current - cell_from_filter = openmc.CellFromFilter(cell1) - cell_to_filter = openmc.CellFilter(cell2) - - cell_to_cell_tallies.append(openmc.Tally(tally_id=tally_index, name=str(cell1)+'-'+str(cell2))) - cell_to_cell_tallies[tally_index].filters = [cell_from_filter, cell_to_filter, energy_filter, polar_filter, azimuthal_filter] - cell_to_cell_tallies[tally_index].scores = ['current'] - tallies_file.append(cell_to_cell_tallies[tally_index]) - tally_index += 1 - - # Cell from + surface filters for partial current - surface_filter = openmc.SurfaceFilter([1]) - cell_to_cell_tallies.append(openmc.Tally(tally_id=tally_index, name=str(cell1)+'-surface1')) - cell_to_cell_tallies[tally_index].filters = [cell_from_filter, surface_filter, energy_filter, polar_filter, azimuthal_filter] - cell_to_cell_tallies[tally_index].scores = ['current'] - tallies_file.append(cell_to_cell_tallies[tally_index]) - tally_index += 1 - - # Surface filter on inner surface, for net current + energy_filter = openmc.EnergyFilter(two_groups) + polar_filter = openmc.PolarFilter([0, np.pi / 4, np.pi]) + azimuthal_filter = openmc.AzimuthalFilter([0, np.pi / 4, np.pi]) surface_filter = openmc.SurfaceFilter([1]) - surf_tally1 = openmc.Tally(tally_id=tally_index, name='net_cylinder') - surf_tally1.filters = [surface_filter, energy_filter, polar_filter, azimuthal_filter] - surf_tally1.scores = ['current'] + cell_from_filter = openmc.CellFromFilter(fuel) + cell_filter = openmc.CellFilter(water) + + # Use Cell to cell filters for partial current + cell_to_cell_tally = openmc.Tally(name=str('fuel_to_water_1')) + cell_to_cell_tally.filters = [cell_from_filter, cell_filter, \ + energy_filter, polar_filter, azimuthal_filter] + cell_to_cell_tally.scores = ['current'] + tallies_file.append(cell_to_cell_tally) + + # Use a Cell from + surface filters for partial current + cell_to_cell_tally = openmc.Tally(name=str('fuel_to_water_2')) + cell_to_cell_tally.filters = [cell_from_filter, surface_filter, \ + energy_filter, polar_filter, azimuthal_filter] + cell_to_cell_tally.scores = ['current'] + tallies_file.append(cell_to_cell_tally) + + # Create partial current tallies from water to fuel + # Filters + cell_from_filter = openmc.CellFromFilter(water) + cell_filter = openmc.CellFilter(fuel) + + # Cell to cell filters for partial current + cell_to_cell_tally = openmc.Tally(name=str('water_to_fuel_1')) + cell_to_cell_tally.filters = [cell_from_filter, cell_filter, \ + energy_filter, polar_filter, azimuthal_filter] + cell_to_cell_tally.scores = ['current'] + tallies_file.append(cell_to_cell_tally) + + # Cell from + surface filters for partial current + cell_to_cell_tally = openmc.Tally(name=str('water_to_fuel_2')) + cell_to_cell_tally.filters = [cell_from_filter, surface_filter, \ + energy_filter, polar_filter, azimuthal_filter] + cell_to_cell_tally.scores = ['current'] + tallies_file.append(cell_to_cell_tally) + + # Create a net current tally on inner surface using a surface filter + surface_filter = openmc.SurfaceFilter([1]) + surf_tally1 = openmc.Tally(name='net_cylinder') + surf_tally1.filters = [surface_filter, energy_filter, polar_filter, \ + azimuthal_filter] + surf_tally1.scores = ['current'] tallies_file.append(surf_tally1) - tally_index += 1 - # Surface filter on left surface, vacuum BC, for net current = leakage + # Create a net current tally on left surface using a surface filter + # This surface has a vacuum boundary condition, so leakage is tallied surface_filter = openmc.SurfaceFilter([2]) - surf_tally2 = openmc.Tally(tally_id=tally_index, name='leakage_left') - surf_tally2.filters = [surface_filter, energy_filter, polar_filter, azimuthal_filter] - surf_tally2.scores = ['current'] + surf_tally2 = openmc.Tally(name='leakage_left') + surf_tally2.filters = [surface_filter, energy_filter, polar_filter, \ + azimuthal_filter] + surf_tally2.scores = ['current'] tallies_file.append(surf_tally2) - tally_index += 1 - # Surface filter on right surface, reflective, current tally doesnt pick up the zero net current though + # Create a net current tally on right surface using a surface filter + # This surface has a reflective boundary condition, but the zero + # net current is not picked up because particles are only tallied once surface_filter = openmc.SurfaceFilter([3]) - surf_tally3 = openmc.Tally(tally_id=tally_index, name='net_right') - surf_tally3.filters = [surface_filter, energy_filter] - surf_tally3.scores = ['current'] + surf_tally3 = openmc.Tally(name='net_right') + surf_tally3.filters = [surface_filter, energy_filter] + surf_tally3.scores = ['current'] + tallies_file.append(surf_tally3) + + surface_filter = openmc.SurfaceFilter([3]) + surf_tally3 = openmc.Tally(name='net_right') + surf_tally3.filters = [surface_filter, energy_filter] + surf_tally3.scores = ['current'] tallies_file.append(surf_tally3) - tally_index += 1 tallies_file.export_to_xml() @@ -153,12 +173,11 @@ class CreateSurfaceTallyTestHarness(PyAPITestHarness): df = df.append(t.get_pandas_dataframe(), ignore_index=True) # Extract the relevant data as a CSV string. - cols = ('d_material', 'd_nuclide', 'd_variable', 'score', 'mean', - 'std. dev.') + cols = ('mean', 'std. dev.') return df.to_csv(None, columns=cols, index=False, float_format='%.7e') return outstr if __name__ == '__main__': - harness = CreateSurfaceTallyTestHarness('statepoint.10.h5') + harness = SurfaceTallyTestHarness('statepoint.10.h5') harness.main() From 1c122733637d93b825eb1c53b947fd5ecb668063 Mon Sep 17 00:00:00 2001 From: guillaume Date: Sun, 6 Aug 2017 14:06:02 -0400 Subject: [PATCH 50/54] updated test_surface_tally results --- tests/test_surface_tally/inputs_true.dat | 50 ++++++----- tests/test_surface_tally/results_true.dat | 104 +++++++++++----------- 2 files changed, 80 insertions(+), 74 deletions(-) diff --git a/tests/test_surface_tally/inputs_true.dat b/tests/test_surface_tally/inputs_true.dat index 213097e8d3..7eeb5102fb 100644 --- a/tests/test_surface_tally/inputs_true.dat +++ b/tests/test_surface_tally/inputs_true.dat @@ -3,11 +3,11 @@ - - - - - + + + + + @@ -44,10 +44,10 @@ - + 10000 - + 10001 @@ -59,7 +59,7 @@ 0.0 0.785398163397 3.14159265359 - + 1 @@ -68,38 +68,42 @@ 10000 - + 2 - + 3 - - 10003 10004 10000 10001 10002 + + 10004 10005 10000 10001 10002 current - - 10003 10005 10000 10001 10002 + + 10004 10003 10000 10001 10002 current - + 10006 10007 10000 10001 10002 current - - 10006 10005 10000 10001 10002 + + 10006 10003 10000 10001 10002 current - - 10005 10000 10001 10002 + + 10003 10000 10001 10002 current - - 10010 10000 10001 10002 + + 10009 10000 10001 10002 current - - 10011 10000 + + 10010 10000 + current + + + 10010 10000 current diff --git a/tests/test_surface_tally/results_true.dat b/tests/test_surface_tally/results_true.dat index df8c1597f4..8e9f6c3900 100644 --- a/tests/test_surface_tally/results_true.dat +++ b/tests/test_surface_tally/results_true.dat @@ -1,51 +1,53 @@ -d_material,d_nuclide,d_variable,score,mean,std. dev. -,,,current,2.4000000e-02,5.6725460e-03 -,,,current,1.0080000e-01,7.0266161e-03 -,,,current,1.2900000e-01,3.1304952e-03 -,,,current,7.2540000e-01,1.4240943e-02 -,,,current,1.8000000e-03,4.1633320e-04 -,,,current,8.3000000e-03,1.4146063e-03 -,,,current,1.5300000e-02,2.9666667e-03 -,,,current,7.9700000e-02,4.3717528e-03 -,,,current,2.4000000e-02,5.6725460e-03 -,,,current,1.0080000e-01,7.0266161e-03 -,,,current,1.2900000e-01,3.1304952e-03 -,,,current,7.2540000e-01,1.4240943e-02 -,,,current,1.8000000e-03,4.1633320e-04 -,,,current,8.3000000e-03,1.4146063e-03 -,,,current,1.5300000e-02,2.9666667e-03 -,,,current,7.9700000e-02,4.3717528e-03 -,,,current,1.5000000e-03,4.0138649e-04 -,,,current,5.7500000e-02,4.1075811e-03 -,,,current,1.0000000e-02,1.3416408e-03 -,,,current,4.5890000e-01,8.6941232e-03 -,,,current,0.0000000e+00,0.0000000e+00 -,,,current,4.6000000e-03,1.0241528e-03 -,,,current,2.0000000e-04,1.3333333e-04 -,,,current,4.6700000e-02,3.5654203e-03 -,,,current,-1.5000000e-03,4.0138649e-04 -,,,current,-5.7500000e-02,4.1075811e-03 -,,,current,-1.0000000e-02,1.3416408e-03 -,,,current,-4.5890000e-01,8.6941232e-03 -,,,current,0.0000000e+00,0.0000000e+00 -,,,current,-4.6000000e-03,1.0241528e-03 -,,,current,-2.0000000e-04,1.3333333e-04 -,,,current,-4.6700000e-02,3.5654203e-03 -,,,current,2.2500000e-02,5.6651959e-03 -,,,current,4.3300000e-02,4.1714639e-03 -,,,current,1.1900000e-01,2.8635642e-03 -,,,current,2.6650000e-01,1.4426249e-02 -,,,current,1.8000000e-03,4.1633320e-04 -,,,current,3.7000000e-03,1.0333333e-03 -,,,current,1.5100000e-02,2.8691849e-03 -,,,current,3.3000000e-02,1.7448336e-03 -,,,current,0.0000000e+00,0.0000000e+00 -,,,current,-5.2600000e-02,4.8355167e-03 -,,,current,0.0000000e+00,0.0000000e+00 -,,,current,-3.9890000e-01,8.1150751e-03 -,,,current,0.0000000e+00,0.0000000e+00 -,,,current,-4.2000000e-03,7.1180522e-04 -,,,current,0.0000000e+00,0.0000000e+00 -,,,current,-3.9000000e-02,4.1686662e-03 -,,,current,-5.9230000e-01,1.2109730e-02 -,,,current,-5.1900000e-02,3.3381632e-03 +mean,std. dev. +2.2000000e-02,1.8915015e-03 +7.0600000e-02,2.4549270e-03 +1.4900000e-01,4.1150132e-03 +6.1190000e-01,9.6014467e-03 +1.1000000e-03,3.7859389e-04 +5.8000000e-03,9.9777530e-04 +1.0300000e-02,1.5205993e-03 +4.6700000e-02,3.0333333e-03 +2.2000000e-02,1.8915015e-03 +7.0600000e-02,2.4549270e-03 +1.4900000e-01,4.1150132e-03 +6.1190000e-01,9.6014467e-03 +1.1000000e-03,3.7859389e-04 +5.8000000e-03,9.9777530e-04 +1.0300000e-02,1.5205993e-03 +4.6700000e-02,3.0333333e-03 +5.4000000e-03,6.6999171e-04 +4.5000000e-02,1.9436506e-03 +3.2300000e-02,1.4067299e-03 +3.8900000e-01,7.4311656e-03 +0.0000000e+00,0.0000000e+00 +2.6000000e-03,4.7609523e-04 +1.0000000e-04,1.0000000e-04 +1.9300000e-02,1.3337499e-03 +-5.4000000e-03,6.6999171e-04 +-4.5000000e-02,1.9436506e-03 +-3.2300000e-02,1.4067299e-03 +-3.8900000e-01,7.4311656e-03 +0.0000000e+00,0.0000000e+00 +-2.6000000e-03,4.7609523e-04 +-1.0000000e-04,1.0000000e-04 +-1.9300000e-02,1.3337499e-03 +1.6600000e-02,1.9675139e-03 +2.5600000e-02,3.2903900e-03 +1.1670000e-01,3.3100856e-03 +2.2290000e-01,5.8924622e-03 +1.1000000e-03,3.7859389e-04 +3.2000000e-03,6.4635731e-04 +1.0200000e-02,1.5114379e-03 +2.7400000e-02,3.0081371e-03 +0.0000000e+00,0.0000000e+00 +-3.5600000e-02,2.2020193e-03 +0.0000000e+00,0.0000000e+00 +-3.6890000e-01,3.4236108e-03 +0.0000000e+00,0.0000000e+00 +-2.8000000e-03,6.7986927e-04 +0.0000000e+00,0.0000000e+00 +-2.2900000e-02,2.4919872e-03 +-7.9490000e-01,1.2723076e-02 +-3.4600000e-02,2.7414514e-03 +-7.9490000e-01,1.2723076e-02 +-3.4600000e-02,2.7414514e-03 From 368842b5df735f7d00afeb30871da05897e8ec84 Mon Sep 17 00:00:00 2001 From: guillaume Date: Sun, 6 Aug 2017 14:06:29 -0400 Subject: [PATCH 51/54] tiny deletion in tally.f90 to avoid a merge --- src/tally.F90 | 3 --- 1 file changed, 3 deletions(-) diff --git a/src/tally.F90 b/src/tally.F90 index 04af19f516..155a31a3ef 100644 --- a/src/tally.F90 +++ b/src/tally.F90 @@ -3271,9 +3271,6 @@ contains ! Reset filter matches flag filter_matches(:) % bins_present = .false. - ! Reset tally map positioning - position = 0 - end subroutine score_surface_tally !=============================================================================== From a766aa33f82f6105d6fa71e51e30349e795c23c1 Mon Sep 17 00:00:00 2001 From: guillaume Date: Sun, 6 Aug 2017 14:19:09 -0400 Subject: [PATCH 52/54] merged with develop then updated surface tally results --- tests/test_surface_tally/inputs_true.dat | 194 +++++++++++------------ 1 file changed, 97 insertions(+), 97 deletions(-) diff --git a/tests/test_surface_tally/inputs_true.dat b/tests/test_surface_tally/inputs_true.dat index 7eeb5102fb..4d18ff1075 100644 --- a/tests/test_surface_tally/inputs_true.dat +++ b/tests/test_surface_tally/inputs_true.dat @@ -1,109 +1,109 @@ - - - - - - - - + + + + + + + + - - - - - - - - - - - - - - - - - - + + + + + + + + + + + + + + + + + + - eigenvalue - 1000 - 10 - 0 - - - -0.62992 -0.62992 -1 0.62992 0.62992 1 - - + eigenvalue + 1000 + 10 + 0 + + + -0.62992 -0.62992 -1 0.62992 0.62992 1 + + - - 10000 - - - 10001 - - - 0.0 4000000.0 20000000.0 - - - 0.0 0.785398163397 3.14159265359 - - - 0.0 0.785398163397 3.14159265359 - - - 1 - - - 10001 - - - 10000 - - - 2 - - - 3 - - - 10004 10005 10000 10001 10002 - current - - - 10004 10003 10000 10001 10002 - current - - - 10006 10007 10000 10001 10002 - current - - - 10006 10003 10000 10001 10002 - current - - - 10003 10000 10001 10002 - current - - - 10009 10000 10001 10002 - current - - - 10010 10000 - current - - - 10010 10000 - current - + + 13 + + + 14 + + + 0.0 4000000.0 20000000.0 + + + 0.0 0.785398163397 3.14159265359 + + + 0.0 0.785398163397 3.14159265359 + + + 1 + + + 14 + + + 13 + + + 2 + + + 3 + + + 5 6 1 2 3 + current + + + 5 4 1 2 3 + current + + + 7 8 1 2 3 + current + + + 7 4 1 2 3 + current + + + 4 1 2 3 + current + + + 10 1 2 3 + current + + + 11 1 + current + + + 11 1 + current + From fa07ddb857ab452a0eda9d4d693efe59f7faca0b Mon Sep 17 00:00:00 2001 From: guillaume Date: Sun, 6 Aug 2017 15:00:38 -0400 Subject: [PATCH 53/54] pleasing python2 --- tests/test_surface_tally/test_surface_tally.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tests/test_surface_tally/test_surface_tally.py b/tests/test_surface_tally/test_surface_tally.py index 5ebe30d587..928b741000 100644 --- a/tests/test_surface_tally/test_surface_tally.py +++ b/tests/test_surface_tally/test_surface_tally.py @@ -87,7 +87,7 @@ class SurfaceTallyTestHarness(PyAPITestHarness): # Create partial current tallies from fuel to water # Filters - two_groups = np.array([0., 4, 20.]) * 1e6 + two_groups = [0., 4e6, 20e6] energy_filter = openmc.EnergyFilter(two_groups) polar_filter = openmc.PolarFilter([0, np.pi / 4, np.pi]) azimuthal_filter = openmc.AzimuthalFilter([0, np.pi / 4, np.pi]) From 56d8f25f36b2cb70d3a698e0dfc35d349e362aa3 Mon Sep 17 00:00:00 2001 From: guillaume Date: Sun, 6 Aug 2017 15:31:02 -0400 Subject: [PATCH 54/54] pleasing python2: .set_element was not giving the same round offs in python2 and python3 --- tests/test_surface_tally/inputs_true.dat | 22 ++--- tests/test_surface_tally/results_true.dat | 92 +++++++++---------- .../test_surface_tally/test_surface_tally.py | 16 ++-- 3 files changed, 62 insertions(+), 68 deletions(-) diff --git a/tests/test_surface_tally/inputs_true.dat b/tests/test_surface_tally/inputs_true.dat index 4d18ff1075..e66d44273a 100644 --- a/tests/test_surface_tally/inputs_true.dat +++ b/tests/test_surface_tally/inputs_true.dat @@ -12,22 +12,16 @@ - - - - - - + + + + - - - - - - - - + + + + diff --git a/tests/test_surface_tally/results_true.dat b/tests/test_surface_tally/results_true.dat index 8e9f6c3900..51e2d7dc0a 100644 --- a/tests/test_surface_tally/results_true.dat +++ b/tests/test_surface_tally/results_true.dat @@ -1,53 +1,53 @@ mean,std. dev. -2.2000000e-02,1.8915015e-03 -7.0600000e-02,2.4549270e-03 -1.4900000e-01,4.1150132e-03 -6.1190000e-01,9.6014467e-03 -1.1000000e-03,3.7859389e-04 -5.8000000e-03,9.9777530e-04 -1.0300000e-02,1.5205993e-03 -4.6700000e-02,3.0333333e-03 -2.2000000e-02,1.8915015e-03 -7.0600000e-02,2.4549270e-03 -1.4900000e-01,4.1150132e-03 -6.1190000e-01,9.6014467e-03 -1.1000000e-03,3.7859389e-04 -5.8000000e-03,9.9777530e-04 -1.0300000e-02,1.5205993e-03 -4.6700000e-02,3.0333333e-03 -5.4000000e-03,6.6999171e-04 -4.5000000e-02,1.9436506e-03 -3.2300000e-02,1.4067299e-03 -3.8900000e-01,7.4311656e-03 +2.1400000e-02,1.1850926e-03 +7.4700000e-02,2.5649345e-03 +1.6090000e-01,4.4433471e-03 +6.4810000e-01,8.6041979e-03 +2.0000000e-03,3.9440532e-04 +4.2000000e-03,8.2731158e-04 +1.1200000e-02,1.0934146e-03 +4.6400000e-02,2.2715633e-03 +2.1400000e-02,1.1850926e-03 +7.4700000e-02,2.5649345e-03 +1.6090000e-01,4.4433471e-03 +6.4810000e-01,8.6041979e-03 +2.0000000e-03,3.9440532e-04 +4.2000000e-03,8.2731158e-04 +1.1200000e-02,1.0934146e-03 +4.6400000e-02,2.2715633e-03 +5.5000000e-03,5.6273143e-04 +4.4300000e-02,1.8502252e-03 +4.5600000e-02,1.2220202e-03 +4.1780000e-01,6.4460151e-03 0.0000000e+00,0.0000000e+00 -2.6000000e-03,4.7609523e-04 -1.0000000e-04,1.0000000e-04 -1.9300000e-02,1.3337499e-03 --5.4000000e-03,6.6999171e-04 --4.5000000e-02,1.9436506e-03 --3.2300000e-02,1.4067299e-03 --3.8900000e-01,7.4311656e-03 +1.4000000e-03,3.3993463e-04 0.0000000e+00,0.0000000e+00 --2.6000000e-03,4.7609523e-04 --1.0000000e-04,1.0000000e-04 --1.9300000e-02,1.3337499e-03 -1.6600000e-02,1.9675139e-03 -2.5600000e-02,3.2903900e-03 -1.1670000e-01,3.3100856e-03 -2.2290000e-01,5.8924622e-03 -1.1000000e-03,3.7859389e-04 -3.2000000e-03,6.4635731e-04 -1.0200000e-02,1.5114379e-03 -2.7400000e-02,3.0081371e-03 +1.6600000e-02,1.0561986e-03 +-5.5000000e-03,5.6273143e-04 +-4.4300000e-02,1.8502252e-03 +-4.5600000e-02,1.2220202e-03 +-4.1780000e-01,6.4460151e-03 0.0000000e+00,0.0000000e+00 --3.5600000e-02,2.2020193e-03 +-1.4000000e-03,3.3993463e-04 0.0000000e+00,0.0000000e+00 --3.6890000e-01,3.4236108e-03 +-1.6600000e-02,1.0561986e-03 +1.5900000e-02,1.1200198e-03 +3.0400000e-02,3.2734623e-03 +1.1530000e-01,3.8327536e-03 +2.3030000e-01,7.2250336e-03 +2.0000000e-03,3.9440532e-04 +2.8000000e-03,7.8598841e-04 +1.1200000e-02,1.0934146e-03 +2.9800000e-02,2.5811281e-03 0.0000000e+00,0.0000000e+00 --2.8000000e-03,6.7986927e-04 +-3.0000000e-02,1.3743685e-03 0.0000000e+00,0.0000000e+00 --2.2900000e-02,2.4919872e-03 --7.9490000e-01,1.2723076e-02 --3.4600000e-02,2.7414514e-03 --7.9490000e-01,1.2723076e-02 --3.4600000e-02,2.7414514e-03 +-3.4810000e-01,5.9711343e-03 +0.0000000e+00,0.0000000e+00 +-3.1000000e-03,6.9041051e-04 +0.0000000e+00,0.0000000e+00 +-2.1900000e-02,2.4241837e-03 +-9.0250000e-01,2.0941320e-02 +-3.5200000e-02,1.2806248e-03 +-9.0250000e-01,2.0941320e-02 +-3.5200000e-02,1.2806248e-03 diff --git a/tests/test_surface_tally/test_surface_tally.py b/tests/test_surface_tally/test_surface_tally.py index 928b741000..ae525a6838 100644 --- a/tests/test_surface_tally/test_surface_tally.py +++ b/tests/test_surface_tally/test_surface_tally.py @@ -13,16 +13,16 @@ class SurfaceTallyTestHarness(PyAPITestHarness): def _build_inputs(self): # Instantiate some Materials and register the appropriate Nuclides uo2 = openmc.Material(name='UO2 fuel at 2.4% wt enrichment') - uo2.set_density('g/cm3', 10.29769) - uo2.add_element('U', 1., enrichment=2.4) - uo2.add_element('O', 2.) + uo2.set_density('g/cc', 10.0) + uo2.add_nuclide('U238', 1.0) + uo2.add_nuclide('U235', 0.02) + uo2.add_nuclide('O16', 2.0) borated_water = openmc.Material(name='Borated water') - borated_water.set_density('g/cm3', 0.740582) - borated_water.add_element('B', 4.0e-5) - borated_water.add_element('H', 5.0e-2) - borated_water.add_element('O', 2.4e-2) - borated_water.add_s_alpha_beta('c_H_in_H2O') + borated_water.set_density('g/cm3', 1) + borated_water.add_nuclide('B10', 10e-5) + borated_water.add_nuclide('H1', 2.0) + borated_water.add_nuclide('O16', 1.0) # Instantiate a Materials collection and export to XML materials_file = openmc.Materials([uo2, borated_water])