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 ! ====================================================================