diff --git a/include/openmc/tallies/tally.h b/include/openmc/tallies/tally.h index 0022fc896..6ffd07be3 100644 --- a/include/openmc/tallies/tally.h +++ b/include/openmc/tallies/tally.h @@ -20,6 +20,9 @@ class Tally { public: Tally() {} + //---------------------------------------------------------------------------- + // Methods for getting and setting filter/stride data. + const std::vector& filters() {return filters_;} int32_t filters(int i) const {return filters_[i];} @@ -30,13 +33,26 @@ public: int32_t n_filter_bins() const {return n_filter_bins_;} + //---------------------------------------------------------------------------- + // Major public data members. + int type_ {TALLY_VOLUME}; //!< volume, surface current //! Event type that contributes to this tally int estimator_ {ESTIMATOR_TRACKLENGTH}; + //! Whether this tally is currently being updated bool active_ {false}; + //! Index of each nuclide to be tallied. -1 indicates total material. + std::vector nuclides_; + + //! True if this tally has a bin for every nuclide in the problem + bool all_nuclides_ {false}; + + //---------------------------------------------------------------------------- + // Miscellaneous public members. + // We need to have quick access to some filters. The following gives indices // for various filters that could be in the tally or C_NONE if they are not // present. @@ -49,6 +65,9 @@ public: int deriv_ {C_NONE}; //!< Index of a TallyDerivative object for diff tallies. private: + //---------------------------------------------------------------------------- + // Private data. + std::vector filters_; //!< Filter indices in global filters array //! Index strides assigned to each filter to support 1D indexing. diff --git a/src/input_xml.F90 b/src/input_xml.F90 index aaab3a460..70d20041b 100644 --- a/src/input_xml.F90 +++ b/src/input_xml.F90 @@ -949,6 +949,8 @@ contains type(XMLNode), allocatable :: node_trigger_list(:) type(DictEntryCI) :: elem type(TallyDerivative), pointer :: deriv + integer, allocatable :: nuclide_bins(:) + logical :: all_nuclides interface subroutine read_tally_derivatives(node_ptr) bind(C) @@ -1193,6 +1195,7 @@ contains ! ======================================================================= ! READ DATA FOR NUCLIDES + all_nuclides = .false. if (check_for_node(node_tal, "nuclides")) then ! Allocate a temporary string array for nuclides and copy values over @@ -1201,27 +1204,23 @@ contains if (trim(sarray(1)) == 'all') then ! Handle special case all - allocate(t % nuclide_bins(n_nuclides + 1)) + allocate(nuclide_bins(n_nuclides + 1)) ! Set bins to 1, 2, 3, ..., n_nuclides, -1 - t % nuclide_bins(1:n_nuclides) = & - (/ (j, j=1, n_nuclides) /) - t % nuclide_bins(n_nuclides + 1) = -1 - - ! Set number of nuclide bins - t % n_nuclide_bins = n_nuclides + 1 + nuclide_bins(1:n_nuclides) = (/ (j, j=1, n_nuclides) /) + nuclide_bins(n_nuclides + 1) = -1 ! Set flag so we can treat this case specially - t % all_nuclides = .true. + all_nuclides = .true. else ! Any other case, e.g. U-235 Pu-239 n_words = node_word_count(node_tal, "nuclides") - allocate(t % nuclide_bins(n_words)) + allocate(nuclide_bins(n_words)) do j = 1, n_words ! Check if total material was specified if (trim(sarray(j)) == 'total') then - t % nuclide_bins(j) = -1 + nuclide_bins(j) = -1 cycle end if @@ -1236,11 +1235,8 @@ contains end if ! Set bin to index in nuclides array - t % nuclide_bins(j) = nuclide_dict % get(word) + nuclide_bins(j) = nuclide_dict % get(word) end do - - ! Set number of nuclide bins - t % n_nuclide_bins = n_words end if ! Deallocate temporary string array @@ -1249,11 +1245,14 @@ contains else ! No were specified -- create only one bin will be added ! for the total material. - allocate(t % nuclide_bins(1)) - t % nuclide_bins(1) = -1 - t % n_nuclide_bins = 1 + allocate(nuclide_bins(1)) + nuclide_bins(1) = -1 end if + ! Set the tally nuclides array + call t % set_nuclide_bins(nuclide_bins, all_nuclides) + deallocate(nuclide_bins) + ! ======================================================================= ! READ DATA FOR SCORES @@ -1290,7 +1289,7 @@ contains select case (trim(score_name)) case ('flux') ! Prohibit user from tallying flux for an individual nuclide - if (.not. (t % n_nuclide_bins == 1 .and. & + if (.not. (t % n_nuclide_bins() == 1 .and. & t % nuclide_bins(1) == -1)) then call fatal_error("Cannot tally flux for an individual nuclide.") end if @@ -1614,8 +1613,8 @@ contains deriv => tally_deriv_c(t % deriv()) if (deriv % variable == DIFF_NUCLIDE_DENSITY & .or. deriv % variable == DIFF_TEMPERATURE) then - if (any(t % nuclide_bins == -1)) then - if (has_energyout) then + do j = 1, t % n_nuclide_bins() + if (has_energyout .and. t % nuclide_bins(j) == -1) then call fatal_error("Error on tally " // trim(to_str(t % id)) & // ": Cannot use a 'nuclide_density' or 'temperature' & &derivative on a tally with an outgoing energy filter and & @@ -1626,7 +1625,7 @@ contains ! (e.g. pertrubing moderator but only tallying fuel), but this ! case would be hard to check for by only reading inputs. end if - end if + end do end if end if diff --git a/src/output.F90 b/src/output.F90 index 248715358..c4e0c3477 100644 --- a/src/output.F90 +++ b/src/output.F90 @@ -625,7 +625,7 @@ contains ! Write results for this filter bin combination score_index = 0 if (t % n_filters() > 0) indent = indent + 2 - do n = 1, t % n_nuclide_bins + do n = 1, t % n_nuclide_bins() ! Write label for nuclide i_nuclide = t % nuclide_bins(n) if (i_nuclide == -1) then diff --git a/src/state_point.F90 b/src/state_point.F90 index 85df5b6df..6de1db773 100644 --- a/src/state_point.F90 +++ b/src/state_point.F90 @@ -277,8 +277,8 @@ contains end if ! Set up nuclide bin array and then write - allocate(str_array(tally % n_nuclide_bins)) - NUCLIDE_LOOP: do j = 1, tally % n_nuclide_bins + allocate(str_array(tally % n_nuclide_bins())) + NUCLIDE_LOOP: do j = 1, tally % n_nuclide_bins() if (tally % nuclide_bins(j) > 0) then if (run_CE) then i_xs = index(nuclides(tally % nuclide_bins(j)) % name, '.') diff --git a/src/tallies/tally.F90 b/src/tallies/tally.F90 index 226491c0e..12d9b9477 100644 --- a/src/tallies/tally.F90 +++ b/src/tallies/tally.F90 @@ -2114,12 +2114,12 @@ contains ! Check for nuclide bins k = 0 - NUCLIDE_LOOP: do while (k < t % n_nuclide_bins) + NUCLIDE_LOOP: do while (k < t % n_nuclide_bins()) ! Increment the index in the list of nuclide bins k = k + 1 - if (t % all_nuclides) then + if (t % all_nuclides()) then ! In the case that the user has requested to tally all nuclides, we ! can take advantage of the fact that we know exactly how nuclide ! bins correspond to nuclide indices. @@ -2259,7 +2259,7 @@ contains ! Nuclide logic ! Check for nuclide bins - NUCLIDE_LOOP: do k = 1, t % n_nuclide_bins + NUCLIDE_LOOP: do k = 1, t % n_nuclide_bins() ! Get index of nuclide in nuclides array i_nuclide = t % nuclide_bins(k) @@ -2623,13 +2623,13 @@ contains ! ====================================================================== ! Nuclide logic - if (t % all_nuclides) then + if (t % all_nuclides()) then if (p % material /= MATERIAL_VOID) then call score_all_nuclides(p, t, flux * filter_weight, filter_index) end if else - NUCLIDE_BIN_LOOP: do k = 1, t % n_nuclide_bins + NUCLIDE_BIN_LOOP: do k = 1, t % n_nuclide_bins() ! Get index of nuclide in nuclides array i_nuclide = t % nuclide_bins(k) @@ -2780,13 +2780,13 @@ contains ! ====================================================================== ! Nuclide logic - if (t % all_nuclides) then + if (t % all_nuclides()) then if (p % material /= MATERIAL_VOID) then call score_all_nuclides(p, t, flux * filter_weight, filter_index) end if else - NUCLIDE_BIN_LOOP: do k = 1, t % n_nuclide_bins + NUCLIDE_BIN_LOOP: do k = 1, t % n_nuclide_bins() ! Get index of nuclide in nuclides array i_nuclide = t % nuclide_bins(k) diff --git a/src/tallies/tally.cpp b/src/tallies/tally.cpp index a93de78b3..81d248ec8 100644 --- a/src/tallies/tally.cpp +++ b/src/tallies/tally.cpp @@ -333,6 +333,20 @@ openmc_tally_set_filters(int32_t index, int n, const int32_t* indices) return 0; } +extern "C" int +openmc_tally_get_nuclides(int32_t index, int** nuclides, int* n) +{ + // Make sure the index fits in the array bounds. + if (index < 1 || index > model::tallies.size()) { + set_errmsg("Index in tallies array is out of bounds."); + return OPENMC_E_OUT_OF_BOUNDS; + } + + //TODO: off-by-one + *n = model::tallies[index-1]->nuclides_.size(); + *nuclides = model::tallies[index-1]->nuclides_.data(); +} + //============================================================================== // Fortran compatibility functions //============================================================================== @@ -391,6 +405,22 @@ extern "C" { int32_t tally_get_n_filter_bins_c(Tally* tally) {return tally->n_filter_bins();} + int tally_get_n_nuclide_bins_c(Tally* tally) + {return tally->nuclides_.size();} + + int tally_get_nuclide_bins_c(Tally* tally, int i) + {return tally->nuclides_[i-1];} + + void + tally_set_nuclide_bins_c(Tally* tally, int n, int bins[], bool all_nuclides) + { + tally->nuclides_.clear(); + tally->nuclides_.assign(bins, bins + n); + tally->all_nuclides_ = all_nuclides; + } + + bool tally_get_all_nuclides_c(Tally* tally) {return tally->all_nuclides_;} + int tally_get_energyin_filter_c(Tally* tally) {return tally->energyin_filter_;} diff --git a/src/tallies/tally_header.F90 b/src/tallies/tally_header.F90 index 2e57d243b..ca90c4087 100644 --- a/src/tallies/tally_header.F90 +++ b/src/tallies/tally_header.F90 @@ -40,6 +40,14 @@ module tally_header integer(C_INT) :: err end function + function openmc_tally_get_nuclides(index, nuclides, n) result(err) bind(C) + import C_INT32_T, C_PTR, C_INT + integer(C_INT32_T), value :: index + type(C_PTR), intent(out) :: nuclides + integer(C_INT), intent(out) :: n + integer(C_INT) :: err + end function openmc_tally_get_nuclides + function active_tallies_size() result(size) bind(C) import C_INT integer(C_INT) :: size @@ -101,11 +109,6 @@ module tally_header real(8) :: volume ! volume of region logical :: depletion_rx = .false. ! has depletion reactions, e.g. (n,2n) - ! Individual nuclides to tally - integer :: n_nuclide_bins = 0 - integer, allocatable :: nuclide_bins(:) - logical :: all_nuclides = .false. - ! Values to score, e.g. flux, absorption, etc. integer :: n_score_bins = 0 integer, allocatable :: score_bins(:) @@ -138,6 +141,10 @@ module tally_header procedure :: filter => tally_get_filter procedure :: stride => tally_get_stride procedure :: n_filter_bins => tally_get_n_filter_bins + procedure :: n_nuclide_bins => tally_get_n_nuclide_bins + procedure :: nuclide_bins => tally_get_nuclide_bins + procedure :: set_nuclide_bins => tally_set_nuclide_bins + procedure :: all_nuclides => tally_get_all_nuclides procedure :: energyin_filter => tally_get_energyin_filter procedure :: energyout_filter => tally_get_energyout_filter procedure :: delayedgroup_filter => tally_get_delayedgroup_filter @@ -289,16 +296,15 @@ contains subroutine tally_allocate_results(this) class(TallyObject), intent(inout) :: this + integer, parameter :: default_nuclide_bins(1) = (/-1/) ! If no nuclides were specified, add a single bin for total material - if (.not. allocated(this % nuclide_bins)) then - allocate(this % nuclide_bins(1)) - this % nuclide_bins(1) = -1 - this % n_nuclide_bins = 1 + if (this % n_nuclide_bins() == 0) then + call this % set_nuclide_bins(default_nuclide_bins) end if ! Set total number of filter and scoring bins - this % total_score_bins = this % n_score_bins * this % n_nuclide_bins + this % total_score_bins = this % n_score_bins * this % n_nuclide_bins() if (allocated(this % results)) then ! If results was already allocated but shape is wrong, then reallocate it @@ -422,6 +428,68 @@ contains n_filter_bins = tally_get_n_filter_bins_c(this % ptr) end function + function tally_get_n_nuclide_bins(this) result(n) + class(TallyObject) :: this + integer(C_INT) :: n + interface + function tally_get_n_nuclide_bins_c(tally) result(n) bind(C) + import C_PTR, C_INT + type(C_PTR), value :: tally + integer(C_INT) :: n + end function + end interface + n = tally_get_n_nuclide_bins_c(this % ptr) + end function + + function tally_get_nuclide_bins(this, i) result(nuclide) + class(TallyObject) :: this + integer(C_INT) :: i, nuclide + interface + function tally_get_nuclide_bins_c(tally, i) result(nuclide) bind(C) + import C_PTR, C_INT + type(C_PTR), value :: tally + integer(C_INT), value :: i + integer(C_INT) :: nuclide + end function + end interface + nuclide = tally_get_nuclide_bins_c(this % ptr, i) + end function + + subroutine tally_set_nuclide_bins(this, bins, all_nuclides) + class(TallyObject) :: this + integer :: bins(:) + logical, optional :: all_nuclides + logical(C_BOOL) :: all_ + interface + subroutine tally_set_nuclide_bins_c(this, n, bins, all_nuclides) bind(C) + import C_PTR, C_INT, C_BOOL + type(C_PTR), value :: this + integer(C_INT), value :: n + integer(C_INT) :: bins(n) + logical(C_BOOL), value :: all_nuclides + end subroutine + end interface + if (present(all_nuclides)) then + all_ = logical(all_nuclides, kind=C_BOOL) + else + all_ = .false._C_BOOL + end if + call tally_set_nuclide_bins_c(this % ptr, size(bins), bins, all_) + end subroutine + + function tally_get_all_nuclides(this) result(all_nuc) + class(TallyObject) :: this + logical(C_BOOL) :: all_nuc + interface + function tally_get_all_nuclides_c(tally) result(all_nuc) bind(C) + import C_PTR, C_BOOL + type(C_PTR), value :: tally + logical(C_BOOL) :: all_nuc + end function + end interface + all_nuc = tally_get_all_nuclides_c(this % ptr) + end function + function tally_get_energyin_filter(this) result(filt) class(TallyObject) :: this integer(C_INT) :: filt @@ -689,31 +757,6 @@ contains end function openmc_tally_get_n_realizations - function openmc_tally_get_nuclides(index, nuclides, n) result(err) bind(C) - ! Return the list of nuclides assigned to a tally - integer(C_INT32_T), value :: index - type(C_PTR), intent(out) :: nuclides - integer(C_INT), intent(out) :: n - integer(C_INT) :: err - - if (index >= 1 .and. index <= size(tallies)) then - associate (t => tallies(index) % obj) - if (allocated(t % nuclide_bins)) then - nuclides = C_LOC(t % nuclide_bins(1)) - n = size(t % nuclide_bins) - err = 0 - else - err = E_ALLOCATE - call set_errmsg("Tally nuclides have not been allocated yet.") - end if - end associate - else - err = E_OUT_OF_BOUNDS - call set_errmsg('Index in tallies array is out of bounds.') - end if - end function openmc_tally_get_nuclides - - function openmc_tally_get_scores(index, scores, n) result(err) bind(C) ! Return the list of nuclides assigned to a tally integer(C_INT32_T), value :: index @@ -856,6 +899,7 @@ contains type(C_PTR), intent(in) :: nuclides(n) integer(C_INT) :: err + integer, allocatable :: bins(:) integer :: i character(C_CHAR), pointer :: string(:) character(len=:, kind=C_CHAR), allocatable :: nuclide_ @@ -863,9 +907,7 @@ contains err = E_UNASSIGNED if (index >= 1 .and. index <= size(tallies)) then associate (t => tallies(index) % obj) - if (allocated(t % nuclide_bins)) deallocate(t % nuclide_bins) - allocate(t % nuclide_bins(n)) - t % n_nuclide_bins = n + allocate(bins(n)) do i = 1, n ! Convert C string to Fortran string @@ -874,10 +916,10 @@ contains select case (nuclide_) case ('total') - t % nuclide_bins(i) = -1 + bins(i) = -1 case default if (nuclide_dict % has(nuclide_)) then - t % nuclide_bins(i) = nuclide_dict % get(nuclide_) + bins(i) = nuclide_dict % get(nuclide_) else err = E_DATA call set_errmsg("Nuclide '" // trim(to_f_string(string)) // & @@ -887,6 +929,8 @@ contains end select end do + call t % set_nuclide_bins(bins) + err = 0 end associate else diff --git a/src/tallies/trigger.F90 b/src/tallies/trigger.F90 index 4c1efcc2e..56ed00aba 100644 --- a/src/tallies/trigger.F90 +++ b/src/tallies/trigger.F90 @@ -183,7 +183,7 @@ contains score_index = trigger % score_index ! Initialize score bin index - NUCLIDE_LOOP: do n = 1, t % n_nuclide_bins + NUCLIDE_LOOP: do n = 1, t % n_nuclide_bins() call get_trigger_uncertainty(std_dev, rel_err, & score_index, filter_index, t)