From ed628b96d2e91743ee614af0fdf1c0cc0eabf0c0 Mon Sep 17 00:00:00 2001 From: Sterling Harper Date: Tue, 5 Feb 2019 01:54:12 -0500 Subject: [PATCH] Move tally scores to C++ --- include/openmc/tallies/tally.h | 10 + src/input_xml.F90 | 281 +------------------------ src/output.F90 | 2 +- src/state_point.F90 | 6 +- src/tallies/tally.F90 | 16 +- src/tallies/tally.cpp | 374 +++++++++++++++++++++++++++++++-- src/tallies/tally_header.F90 | 250 +++++----------------- src/tallies/trigger.cpp | 7 +- 8 files changed, 432 insertions(+), 514 deletions(-) diff --git a/include/openmc/tallies/tally.h b/include/openmc/tallies/tally.h index 6c27d9655a..b9b97f64e2 100644 --- a/include/openmc/tallies/tally.h +++ b/include/openmc/tallies/tally.h @@ -21,6 +21,10 @@ class Tally { public: Tally() {} + void set_scores(pugi::xml_node node); + + void set_scores(std::vector scores); + //---------------------------------------------------------------------------- // Methods for getting and setting filter/stride data. @@ -37,6 +41,8 @@ public: //---------------------------------------------------------------------------- // Other methods. + void init_scores(pugi::xml_node node); + void init_triggers(pugi::xml_node node, int i_tally); //---------------------------------------------------------------------------- @@ -52,6 +58,8 @@ public: //! Whether this tally is currently being updated bool active_ {false}; + std::vector scores_; //!< Filter integrands (e.g. flux, fission) + //! Index of each nuclide to be tallied. -1 indicates total material. std::vector nuclides_; @@ -67,6 +75,8 @@ public: int energyout_filter_ {C_NONE}; int delayedgroup_filter_ {C_NONE}; + bool depletion_rx_ {false}; //!< Has depletion reactions (e.g. (n,2n)) + std::vector triggers_; int deriv_ {C_NONE}; //!< Index of a TallyDerivative object for diff tallies. diff --git a/src/input_xml.F90 b/src/input_xml.F90 index a9192ba3e4..e1aa1b5033 100644 --- a/src/input_xml.F90 +++ b/src/input_xml.F90 @@ -943,6 +943,12 @@ contains logical :: all_nuclides interface + subroutine tally_init_scores(tally_ptr, xml_node) bind(C) + import C_PTR + type(C_PTR), value :: tally_ptr + type(C_PTR) :: xml_node + end subroutine + subroutine tally_init_triggers(tally_ptr, i_tally, xml_node) bind(C) import C_PTR, C_INT type(C_PTR), value :: tally_ptr @@ -1255,281 +1261,10 @@ contains ! ======================================================================= ! READ DATA FOR SCORES + call tally_init_scores(t % ptr, node_tal % ptr) + if (check_for_node(node_tal, "scores")) then n_words = node_word_count(node_tal, "scores") - allocate(sarray(n_words)) - call get_node_array(node_tal, "scores", sarray) - n_scores = n_words - - ! Allocate score storage accordingly - allocate(t % score_bins(n_scores)) - - ! Check the validity of the scores and their filters - do j = 1, n_scores - score_name = sarray(j) - - ! Check if delayed group filter is used with any score besides - ! delayed-nu-fission or decay-rate - if ((score_name /= 'delayed-nu-fission' .and. & - score_name /= 'decay-rate') .and. has_delayedgroup) then - call fatal_error("Cannot tally " // trim(score_name) // " with a & - &delayedgroup filter.") - end if - - select case (trim(score_name)) - case ('flux') - ! Prohibit user from tallying flux for an individual nuclide - 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 - - t % score_bins(j) = SCORE_FLUX - if (has_energyout) then - call fatal_error("Cannot tally flux with an outgoing energy & - &filter.") - end if - - case ('total', '(n,total)') - t % score_bins(j) = SCORE_TOTAL - if (has_energyout) then - call fatal_error("Cannot tally total reaction rate with an & - &outgoing energy filter.") - end if - - case ('scatter') - t % score_bins(j) = SCORE_SCATTER - if (has_energyout .or. has_legendre) then - ! Set tally estimator to analog - call t % set_estimator(ESTIMATOR_ANALOG) - end if - - case ('nu-scatter') - t % score_bins(j) = SCORE_NU_SCATTER - - ! Set tally estimator to analog for CE mode - ! (MG mode has all data available without a collision being - ! necessary) - if (run_CE) then - call t % set_estimator(ESTIMATOR_ANALOG) - else - if (has_energyout .or. has_legendre) then - ! Set tally estimator to analog - call t % set_estimator(ESTIMATOR_ANALOG) - end if - end if - - case ('n2n', '(n,2n)') - t % score_bins(j) = N_2N - t % depletion_rx = .true. - - case ('n3n', '(n,3n)') - t % score_bins(j) = N_3N - t % depletion_rx = .true. - - case ('n4n', '(n,4n)') - t % score_bins(j) = N_4N - t % depletion_rx = .true. - - case ('absorption') - t % score_bins(j) = SCORE_ABSORPTION - if (has_energyout) then - call fatal_error("Cannot tally absorption rate with an outgoing & - &energy filter.") - end if - case ('fission', '18') - t % score_bins(j) = SCORE_FISSION - if (has_energyout) then - call fatal_error("Cannot tally fission rate with an outgoing & - &energy filter.") - end if - case ('nu-fission') - t % score_bins(j) = SCORE_NU_FISSION - if (has_energyout) then - ! Set tally estimator to analog - call t % set_estimator(ESTIMATOR_ANALOG) - end if - case ('decay-rate') - t % score_bins(j) = SCORE_DECAY_RATE - case ('delayed-nu-fission') - t % score_bins(j) = SCORE_DELAYED_NU_FISSION - if (has_energyout) then - ! Set tally estimator to analog - call t % set_estimator(ESTIMATOR_ANALOG) - end if - case ('prompt-nu-fission') - t % score_bins(j) = SCORE_PROMPT_NU_FISSION - if (has_energyout) then - ! Set tally estimator to analog - call t % set_estimator(ESTIMATOR_ANALOG) - end if - case ('kappa-fission') - t % score_bins(j) = SCORE_KAPPA_FISSION - case ('inverse-velocity') - t % score_bins(j) = SCORE_INVERSE_VELOCITY - case ('fission-q-prompt') - t % score_bins(j) = SCORE_FISS_Q_PROMPT - case ('fission-q-recoverable') - t % score_bins(j) = SCORE_FISS_Q_RECOV - case ('current') - - ! Check which type of current is desired: mesh currents or - ! surface currents - if (has_surface .or. has_cell .or. has_cellfrom) then - - ! Check to make sure that mesh surface currents are not desired as well - if (has_meshsurface) then - call fatal_error("Cannot tally mesh surface currents & - &in the same tally as normal surface currents") - end if - - call t % set_type(TALLY_SURFACE) - t % score_bins(j) = SCORE_CURRENT - - else if (has_meshsurface) then - t % score_bins(j) = SCORE_CURRENT - call t % set_type(TALLY_MESH_SURFACE) - - ! 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 - else - call fatal_error("Cannot tally currents without surface & - &type filters") - end if - - case ('events') - t % score_bins(j) = SCORE_EVENTS - case ('elastic', '(n,elastic)') - t % score_bins(j) = ELASTIC - case ('(n,2nd)') - t % score_bins(j) = N_2ND - case ('(n,na)') - t % score_bins(j) = N_2NA - case ('(n,n3a)') - t % score_bins(j) = N_N3A - case ('(n,2na)') - t % score_bins(j) = N_2NA - case ('(n,3na)') - t % score_bins(j) = N_3NA - case ('(n,np)') - t % score_bins(j) = N_NP - case ('(n,n2a)') - t % score_bins(j) = N_N2A - case ('(n,2n2a)') - t % score_bins(j) = N_2N2A - case ('(n,nd)') - t % score_bins(j) = N_ND - case ('(n,nt)') - t % score_bins(j) = N_NT - case ('(n,nHe-3)') - t % score_bins(j) = N_N3HE - case ('(n,nd2a)') - t % score_bins(j) = N_ND2A - case ('(n,nt2a)') - t % score_bins(j) = N_NT2A - case ('(n,3nf)') - t % score_bins(j) = N_3NF - case ('(n,2np)') - t % score_bins(j) = N_2NP - case ('(n,3np)') - t % score_bins(j) = N_3NP - case ('(n,n2p)') - t % score_bins(j) = N_N2P - case ('(n,npa)') - t % score_bins(j) = N_NPA - case ('(n,n1)') - t % score_bins(j) = N_N1 - case ('(n,nc)') - t % score_bins(j) = N_NC - case ('(n,gamma)') - t % score_bins(j) = N_GAMMA - t % depletion_rx = .true. - case ('(n,p)') - t % score_bins(j) = N_P - t % depletion_rx = .true. - case ('(n,d)') - t % score_bins(j) = N_D - case ('(n,t)') - t % score_bins(j) = N_T - case ('(n,3He)') - t % score_bins(j) = N_3HE - case ('(n,a)') - t % score_bins(j) = N_A - t % depletion_rx = .true. - case ('(n,2a)') - t % score_bins(j) = N_2A - case ('(n,3a)') - t % score_bins(j) = N_3A - case ('(n,2p)') - t % score_bins(j) = N_2P - case ('(n,pa)') - t % score_bins(j) = N_PA - case ('(n,t2a)') - t % score_bins(j) = N_T2A - case ('(n,d2a)') - t % score_bins(j) = N_D2A - case ('(n,pd)') - t % score_bins(j) = N_PD - case ('(n,pt)') - t % score_bins(j) = N_PT - case ('(n,da)') - t % score_bins(j) = N_DA - - case default - ! First look for deprecated scores - if (starts_with(trim(score_name), 'scatter-') .or. & - starts_with(trim(score_name), 'nu-scatter-') .or. & - starts_with(trim(score_name), 'total-y') .or. & - starts_with(trim(score_name), 'flux-y')) then - call fatal_error(trim(score_name) // " is no longer available.") - end if - - ! Assume that user has specified an MT number - MT = int(str_to_int(score_name)) - - if (MT /= ERROR_INT) then - ! Specified score was an integer - if (MT > 1) then - t % score_bins(j) = MT - else - call fatal_error("Invalid MT on : " // trim(score_name)) - end if - - else - ! Specified score was not an integer - call fatal_error("Unknown scoring function: " // trim(score_name)) - end if - - end select - - ! Do a check at the end (instead of for every case) to make sure - ! the tallies are compatible with MG mode where we have less detailed - ! nuclear data - if (.not. run_CE .and. t % score_bins(j) > 0) then - call fatal_error("Cannot tally " // trim(score_name) // & - " reaction rate in multi-group mode") - end if - end do - - t % n_score_bins = n_scores - - ! Deallocate temporary string array of scores - deallocate(sarray) - - ! Check that no duplicate scores exist - do j = 1, n_scores - 1 - do k = j + 1, n_scores - if (t % score_bins(j) == t % score_bins(k)) then - call fatal_error("Duplicate score of type '" // trim(& - reaction_name(t % score_bins(j))) // "' found in tally " & - // trim(to_str(t % id()))) - end if - end do - end do ! Check if tally is compatible with particle type if (photon_transport) then diff --git a/src/output.F90 b/src/output.F90 index 7bbf2094db..0a3781599e 100644 --- a/src/output.F90 +++ b/src/output.F90 @@ -641,7 +641,7 @@ contains end if indent = indent + 2 - do k = 1, t % n_score_bins + do k = 1, t % n_score_bins() score_index = score_index + 1 associate(r => t % results(RESULT_SUM:RESULT_SUM_SQ, :, :)) diff --git a/src/state_point.F90 b/src/state_point.F90 index fca7954ab0..72c2e090b6 100644 --- a/src/state_point.F90 +++ b/src/state_point.F90 @@ -312,9 +312,9 @@ contains end if ! Write scores. - call write_dataset(tally_group, "n_score_bins", tally % n_score_bins) - allocate(str_array(size(tally % score_bins))) - do j = 1, size(tally % score_bins) + call write_dataset(tally_group, "n_score_bins", tally % n_score_bins()) + allocate(str_array(tally % n_score_bins())) + do j = 1, tally % n_score_bins() str_array(j) = reaction_name(tally % score_bins(j)) end do call write_dataset(tally_group, "score_bins", str_array) diff --git a/src/tallies/tally.F90 b/src/tallies/tally.F90 index 05ee7ac92f..9d45eae87f 100644 --- a/src/tallies/tally.F90 +++ b/src/tallies/tally.F90 @@ -142,7 +142,7 @@ contains ! Pre-collision energy of particle E = p % last_E - SCORE_LOOP: do i = 1, t % n_score_bins + SCORE_LOOP: do i = 1, t % n_score_bins() ! determine what type of score bin score_bin = t % score_bins(i) @@ -1308,7 +1308,7 @@ contains end if i = 0 - SCORE_LOOP: do q = 1, t % n_score_bins + SCORE_LOOP: do q = 1, t % n_score_bins() i = i + 1 ! determine what type of score bin @@ -2064,7 +2064,7 @@ contains atom_density = mat % atom_density(i) ! Determine score for each bin - call score_general(p, i_tally, (i_nuclide-1)*t % n_score_bins, filter_index, & + call score_general(p, i_tally, (i_nuclide-1)*t % n_score_bins(), filter_index, & i_nuclide, atom_density, flux) end do NUCLIDE_LOOP @@ -2076,7 +2076,7 @@ contains atom_density = ZERO ! Determine score for each bin - call score_general(p, i_tally, n_nuclides*t % n_score_bins, filter_index, & + call score_general(p, i_tally, n_nuclides*t % n_score_bins(), filter_index, & i_nuclide, atom_density, flux) end associate @@ -3123,7 +3123,7 @@ contains err = openmc_tally_get_active(i, active) if (active) then ! Check if tally contains depletion reactions and if so, set flag - if (t % depletion_rx) need_depletion_rx = .true. + if (t % depletion_rx()) need_depletion_rx = .true. end if end associate end do @@ -3196,10 +3196,4 @@ contains dens = materials(i_material) % atom_density(i) end function - function tally_get_n_score_bins(i_tally) result(n) bind(C) - integer(C_INT), value :: i_tally - integer(C_INT) :: n - n = tallies(i_tally) % obj % n_score_bins - end function - end module tally diff --git a/src/tallies/tally.cpp b/src/tallies/tally.cpp index 2e8bccc3cb..1199f6ba62 100644 --- a/src/tallies/tally.cpp +++ b/src/tallies/tally.cpp @@ -8,10 +8,14 @@ #include "openmc/settings.h" #include "openmc/tallies/derivative.h" #include "openmc/tallies/filter.h" -#include "openmc/tallies/filter_energy.h" +#include "openmc/tallies/filter_cell.h" +#include "openmc/tallies/filter_cellfrom.h" #include "openmc/tallies/filter_delayedgroup.h" -#include "openmc/tallies/filter_surface.h" +#include "openmc/tallies/filter_energy.h" +#include "openmc/tallies/filter_legendre.h" #include "openmc/tallies/filter_mesh.h" +#include "openmc/tallies/filter_meshsurface.h" +#include "openmc/tallies/filter_surface.h" #include "openmc/xml_interface.h" #include "xtensor/xadapt.hpp" @@ -44,8 +48,6 @@ extern "C" int material_nuclide_index(int i_material, int i_nuclide); extern "C" double material_atom_density(int i_material, int i); -extern "C" int tally_get_n_score_bins(int i_tally); - //============================================================================== // Global variable definitions //============================================================================== @@ -189,6 +191,162 @@ private: const Tally& tally_; }; +int +score_str_to_int(std::string score_str) +{ + if (score_str == "flux") + return SCORE_FLUX; + + if (score_str == "total" || score_str == "(n,total)") + return SCORE_TOTAL; + + if (score_str == "scatter") + return SCORE_SCATTER; + + if (score_str == "nu-scatter") + return SCORE_NU_SCATTER; + + if (score_str == "absorption") + return SCORE_ABSORPTION; + + if (score_str == "fission" || score_str == "18") + return SCORE_FISSION; + + if (score_str == "nu-fission") + return SCORE_NU_FISSION; + + if (score_str == "decay-rate") + return SCORE_DECAY_RATE; + + if (score_str == "delayed-nu-fission") + return SCORE_DELAYED_NU_FISSION; + + if (score_str == "prompt-nu-fission") + return SCORE_PROMPT_NU_FISSION; + + if (score_str == "kappa-fission") + return SCORE_KAPPA_FISSION; + + if (score_str == "inverse-velocity") + return SCORE_INVERSE_VELOCITY; + + if (score_str == "fission-q-prompt") + return SCORE_FISS_Q_PROMPT; + + if (score_str == "fission-q-recoverable") + return SCORE_FISS_Q_RECOV; + + if (score_str == "current") + return SCORE_CURRENT; + + if (score_str == "events") + return SCORE_EVENTS; + + if (score_str == "elastic" || score_str == "(n,elastic)") + return ELASTIC; + + if (score_str == "n2n" || score_str == "(n,2n)") + return N_2N; + + if (score_str == "n3n" || score_str == "(n,3n)") + return N_3N; + + if (score_str == "n4n" || score_str == "(n,4n)") + return N_4N; + + if (score_str == "(n,2nd)") + return N_2ND; + if (score_str == "(n,na)") + return N_2NA; + if (score_str == "(n,n3a)") + return N_N3A; + if (score_str == "(n,2na)") + return N_2NA; + if (score_str == "(n,3na)") + return N_3NA; + if (score_str == "(n,np)") + return N_NP; + if (score_str == "(n,n2a)") + return N_N2A; + if (score_str == "(n,2n2a)") + return N_2N2A; + if (score_str == "(n,nd)") + return N_ND; + if (score_str == "(n,nt)") + return N_NT; + if (score_str == "(n,nHe-3)") + return N_N3HE; + if (score_str == "(n,nd2a)") + return N_ND2A; + if (score_str == "(n,nt2a)") + return N_NT2A; + if (score_str == "(n,3nf)") + return N_3NF; + if (score_str == "(n,2np)") + return N_2NP; + if (score_str == "(n,3np)") + return N_3NP; + if (score_str == "(n,n2p)") + return N_N2P; + if (score_str == "(n,npa)") + return N_NPA; + if (score_str == "(n,n1)") + return N_N1; + if (score_str == "(n,nc)") + return N_NC; + if (score_str == "(n,gamma)") + return N_GAMMA; + if (score_str == "(n,p)") + return N_P; + if (score_str == "(n,d)") + return N_D; + if (score_str == "(n,t)") + return N_T; + if (score_str == "(n,3He)") + return N_3HE; + if (score_str == "(n,a)") + return N_A; + if (score_str == "(n,2a)") + return N_2A; + if (score_str == "(n,3a)") + return N_3A; + if (score_str == "(n,2p)") + return N_2P; + if (score_str == "(n,pa)") + return N_PA; + if (score_str == "(n,t2a)") + return N_T2A; + if (score_str == "(n,d2a)") + return N_D2A; + if (score_str == "(n,pd)") + return N_PD; + if (score_str == "(n,pt)") + return N_PT; + if (score_str == "(n,da)") + return N_DA; + + // So far we have not identified this score string. Check to see if it is a + // deprecated score. + if (score_str.rfind("scatter-", 0) == 0 + || score_str.rfind("nu-scatter-", 0) == 0 + || score_str.rfind("total-y", 0) == 0 + || score_str.rfind("flux-y", 0) == 0) + fatal_error(score_str + " is no longer an available score"); + + + // Assume the given string is a reaction MT number. Make sure it's a natural + // number then return. + int MT; + try { + MT = std::stoi(score_str); + } catch (const std::invalid_argument& ex) { + throw std::invalid_argument("Invalid tally score \"" + score_str + "\""); + } + if (MT < 1) + throw std::invalid_argument("Invalid tally score \"" + score_str + "\""); + return MT; +} + //============================================================================== // Tally object implementation //============================================================================== @@ -233,6 +391,147 @@ Tally::set_filters(const int32_t filter_indices[], int n) n_filter_bins_ = stride; } +void +Tally::set_scores(pugi::xml_node node) +{ + if (!check_for_node(node, "scores")) + fatal_error("No scores specified on tally " + std::to_string(id_)); + + auto scores = get_node_array(node, "scores"); + set_scores(scores); +} + +void +Tally::set_scores(std::vector scores) +{ + // Reset state and prepare for the new scores. + scores_.clear(); + depletion_rx_ = false; + scores_.reserve(scores.size()); + + // Check for the presence of certain restrictive filters. + bool energyout_present = energyout_filter_ != C_NONE; + //bool delayedgroup_present = false; + bool legendre_present = false; + bool cell_present = false; + bool cellfrom_present = false; + bool surface_present = false; + bool meshsurface_present = false; + for (auto i_filt : filters_) { + //TODO: off-by-one + const auto* filt {model::tally_filters[i_filt-1].get()}; + if (dynamic_cast(filt)) { + legendre_present = true; + } else if (dynamic_cast(filt)) { + cellfrom_present = true; + } else if (dynamic_cast(filt)) { + cell_present = true; + } else if (dynamic_cast(filt)) { + surface_present = true; + } else if (dynamic_cast(filt)) { + meshsurface_present = true; + } + } + + // Iterate over the given scores. + for (auto score_str : scores) { + // Make sure a delayed group filter wasn't used with an incompatible score. + bool has_delayedgroup = delayedgroup_filter_ != C_NONE; + if (delayedgroup_filter_ != C_NONE) { + if (score_str != "delayed-nu-fission" && score_str != "decay-rate") + fatal_error("Cannot tally " + score_str + "with a delayedgroup filter"); + } + + auto score = score_str_to_int(score_str); + + switch (score) { + case SCORE_FLUX: + if (!nuclides_.empty()) + if (!(nuclides_.size() == 1 && nuclides_[0] == -1)) + fatal_error("Cannot tally flux for an individual nuclide."); + if (energyout_present) + fatal_error("Cannot tally flux with an outgoing energy filter."); + break; + + case SCORE_TOTAL: + case SCORE_ABSORPTION: + case SCORE_FISSION: + if (energyout_present) + fatal_error("Cannot tally " + score_str + " reaction rate with an " + "outgoing energy filter"); + break; + + case SCORE_SCATTER: + if (legendre_present) + estimator_ = ESTIMATOR_ANALOG; + case SCORE_NU_FISSION: + case SCORE_DELAYED_NU_FISSION: + case SCORE_PROMPT_NU_FISSION: + if (energyout_present) + estimator_ = ESTIMATOR_ANALOG; + break; + + case SCORE_NU_SCATTER: + if (settings::run_CE) { + estimator_ = ESTIMATOR_ANALOG; + } else { + if (energyout_present || legendre_present) + estimator_ = ESTIMATOR_ANALOG; + } + break; + + case N_2N: + case N_3N: + case N_4N: + case N_GAMMA: + case N_P: + case N_A: + depletion_rx_ = true; + break; + + case SCORE_CURRENT: + // Check which type of current is desired: mesh or surface currents. + if (surface_present || cell_present || cellfrom_present) { + if (meshsurface_present) + fatal_error("Cannot tally mesh surface currents in the same tally as " + "normal surface currents"); + type_ = TALLY_SURFACE; + } else if (meshsurface_present) { + type_ = TALLY_MESH_SURFACE; + } else { + fatal_error("Cannot tally currents without surface type filters"); + } + break; + } + + scores_.push_back(score); + } + + // Make sure that no duplicate scores exist. + for (auto it1 = scores_.begin(); it1 != scores_.end(); ++it1) { + for (auto it2 = it1 + 1; it2 != scores_.end(); ++it2) { + if (*it1 == *it2) + fatal_error("Duplicate score of type \"" + reaction_name(*it1) + + "\" found in tally " + std::to_string(id_)); + } + } + + // Make sure all scores are compatible with multigroup mode. + if (!settings::run_CE) { + for (auto sc : scores_) + if (sc > 0) + fatal_error("Cannot tally " + reaction_name(sc) + " reaction rate " + "in multi-group mode"); + } + + // Make sure current scores are not mixed in with volumetric scores. + if (type_ == TALLY_SURFACE || type_ == TALLY_MESH_SURFACE) { + if (scores_.size() != 1) + fatal_error("Cannot tally other scores in the same tally as surface " + "currents"); + } +} + void Tally::init_triggers(pugi::xml_node node, int i_tally) { @@ -347,7 +646,6 @@ score_analog_tally_ce(Particle* p) for (auto i_tally : model::active_analog_tallies) { //TODO: off-by-one const Tally& tally {*model::tallies[i_tally-1]}; - auto n_score_bins = tally_get_n_score_bins(i_tally); // Initialize an iterator over valid filter bin combinations. If there are // no valid combinations, use a continue statement to ensure we skip the @@ -371,7 +669,7 @@ score_analog_tally_ce(Particle* p) // and flux arguments for score_general are not used for analog // tallies. if (i_nuclide == p->event_nuclide || i_nuclide == -1) - score_general_ce(p, i_tally, i*n_score_bins, filter_index, + score_general_ce(p, i_tally, i*tally.scores_.size(), filter_index, -1, -1., filter_weight); } @@ -380,12 +678,12 @@ score_analog_tally_ce(Particle* p) // can take advantage of the fact that we know exactly how nuclide // bins correspond to nuclide indices. First, tally the nuclide. auto i = p->event_nuclide; - score_general_ce(p, i_tally, i*n_score_bins, filter_index, + score_general_ce(p, i_tally, i*tally.scores_.size(), filter_index, -1, -1., filter_weight); // Now tally the total material. i = tally.nuclides_.size(); - score_general_ce(p, i_tally, i*n_score_bins, filter_index, + score_general_ce(p, i_tally, i*tally.scores_.size(), filter_index, -1, -1., filter_weight); } } @@ -412,7 +710,6 @@ score_analog_tally_mg(Particle* p) for (auto i_tally : model::active_analog_tallies) { //TODO: off-by-one const Tally& tally {*model::tallies[i_tally-1]}; - auto n_score_bins = tally_get_n_score_bins(i_tally); // Initialize an iterator over valid filter bin combinations. If there are // no valid combinations, use a continue statement to ensure we skip the @@ -437,7 +734,7 @@ score_analog_tally_mg(Particle* p) atom_density = material_atom_density(p->material, j); } - score_general_mg(p, i_tally, i*n_score_bins, filter_index, + score_general_mg(p, i_tally, i*tally.scores_.size(), filter_index, i_nuclide, atom_density, filter_weight); } } @@ -469,7 +766,6 @@ score_tracklength_tally(Particle* p, double distance) for (auto i_tally : model::active_tracklength_tallies) { //TODO: off-by-one const Tally& tally {*model::tallies[i_tally-1]}; - auto n_score_bins = tally_get_n_score_bins(i_tally); // Initialize an iterator over valid filter bin combinations. If there are // no valid combinations, use a continue statement to ensure we skip the @@ -503,10 +799,10 @@ score_tracklength_tally(Particle* p, double distance) //TODO: consider replacing this "if" with pointers or templates if (settings::run_CE) { - score_general_ce(p, i_tally, i*n_score_bins, filter_index, + score_general_ce(p, i_tally, i*tally.scores_.size(), filter_index, i_nuclide, atom_density, flux*filter_weight); } else { - score_general_mg(p, i_tally, i*n_score_bins, filter_index, + score_general_mg(p, i_tally, i*tally.scores_.size(), filter_index, i_nuclide, atom_density, flux*filter_weight); } } @@ -548,7 +844,6 @@ score_collision_tally(Particle* p) for (auto i_tally : model::active_collision_tallies) { //TODO: off-by-one const Tally& tally {*model::tallies[i_tally-1]}; - auto n_score_bins = tally_get_n_score_bins(i_tally); // Initialize an iterator over valid filter bin combinations. If there are // no valid combinations, use a continue statement to ensure we skip the @@ -579,10 +874,10 @@ score_collision_tally(Particle* p) //TODO: consider replacing this "if" with pointers or templates if (settings::run_CE) { - score_general_ce(p, i_tally, i*n_score_bins, filter_index, + score_general_ce(p, i_tally, i*tally.scores_.size(), filter_index, i_nuclide, atom_density, flux*filter_weight); } else { - score_general_mg(p, i_tally, i*n_score_bins, filter_index, + score_general_mg(p, i_tally, i*tally.scores_.size(), filter_index, i_nuclide, atom_density, flux*filter_weight); } } @@ -612,7 +907,6 @@ score_surface_tally_inner(Particle* p, const std::vector& tallies) for (auto i_tally : tallies) { //TODO: off-by-one const Tally& tally {*model::tallies[i_tally-1]}; - auto n_score_bins = tally_get_n_score_bins(i_tally); auto results = tally_results(i_tally); // Initialize an iterator over valid filter bin combinations. If there are @@ -631,7 +925,8 @@ score_surface_tally_inner(Particle* p, const std::vector& tallies) // There is only one score type for current tallies so there is no need // for a further scoring function. double score = flux * filter_weight; - for (auto score_index = 1; score_index < n_score_bins+1; ++score_index) { + for (auto score_index = 1; score_index < tally.scores_.size()+1; + ++score_index) { //TODO: off-by-one #pragma omp atomic results(filter_index-1, score_index-1, RESULT_VALUE) += score; @@ -853,6 +1148,40 @@ openmc_tally_set_active(int32_t index, bool active) return 0; } +extern "C" int +openmc_tally_get_scores(int32_t index, int** scores, int* n) +{ + 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 + *scores = model::tallies[index-1]->scores_.data(); + *n = model::tallies[index-1]->scores_.size(); + return 0; +} + +extern "C" int +openmc_tally_set_scores(int32_t index, int n, const char** scores) +{ + if (index < 1 || index > model::tallies.size()) { + set_errmsg("Index in tallies array is out of bounds."); + return OPENMC_E_OUT_OF_BOUNDS; + } + + std::vector scores_str(scores, scores+n); + try { + //TODO: off-by-one + model::tallies[index-1]->set_scores(scores_str); + } catch (const std::invalid_argument& ex) { + set_errmsg(ex.what()); + return OPENMC_E_INVALID_ARGUMENT; + } + + return 0; +} + extern "C" int openmc_tally_get_filters(int32_t index, const int32_t** indices, int* n) { @@ -951,6 +1280,12 @@ extern "C" { void tally_set_estimator_c(Tally* tally, int e) {tally->estimator_ = e;} + bool tally_get_depletion_rx_c(Tally* tally) {return tally->depletion_rx_;} + + int tally_get_n_scores_c(Tally* tally) {return tally->scores_.size();} + + int tally_get_score_c(Tally* tally, int i) {return tally->scores_[i];} + void tally_set_filters_c(Tally* tally, int n, int32_t filter_indices[]) {tally->set_filters(filter_indices, n);} @@ -985,6 +1320,9 @@ extern "C" { int tally_get_delayedgroup_filter_c(Tally* tally) {return tally->delayedgroup_filter_;} + void tally_init_scores(Tally* tally, pugi::xml_node* node) + {tally->set_scores(*node);} + void tally_init_triggers(Tally* tally, int i_tally, pugi::xml_node* node) {tally->init_triggers(*node, i_tally);} diff --git a/src/tallies/tally_header.F90 b/src/tallies/tally_header.F90 index 275aa0df32..60e339e2b7 100644 --- a/src/tallies/tally_header.F90 +++ b/src/tallies/tally_header.F90 @@ -97,11 +97,10 @@ module tally_header character(len=104) :: name = "" ! user-defined name real(8) :: volume ! volume of region - logical :: depletion_rx = .false. ! has depletion reactions, e.g. (n,2n) ! Values to score, e.g. flux, absorption, etc. - integer :: n_score_bins = 0 - integer, allocatable :: score_bins(:) + !integer :: n_score_bins = 0 + !integer, allocatable :: score_bins_(:) ! Results for each bin -- the first dimension of the array is for scores ! (e.g. flux, total reaction rate, fission reaction rate, etc.) and the @@ -125,6 +124,9 @@ module tally_header procedure :: set_type => tally_set_type procedure :: estimator => tally_get_estimator procedure :: set_estimator => tally_set_estimator + procedure :: depletion_rx => tally_get_depletion_rx + procedure :: n_score_bins => tally_get_n_score_bins + procedure :: score_bins => tally_get_score_bin procedure :: n_filters => tally_get_n_filters procedure :: filter => tally_get_filter procedure :: stride => tally_get_stride @@ -285,7 +287,7 @@ contains 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 @@ -379,6 +381,47 @@ contains call tally_set_estimator_c(this % ptr, e) end subroutine + function tally_get_depletion_rx(this) result(drx) + class(TallyObject) :: this + logical(C_BOOL) :: drx + interface + function tally_get_depletion_rx_c(tally) result(drx) bind(C) + import C_PTR, C_BOOL + type(C_PTR), value :: tally + logical(C_BOOl) :: drx + end function + end interface + drx = tally_get_depletion_rx_c(this % ptr) + end function + + function tally_get_n_score_bins(this) result(n) + class(TallyObject) :: this + integer(C_INT) :: n + interface + function tally_get_n_scores_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_scores_c(this % ptr) + end function + + function tally_get_score_bin(this, i) result(filt) + class(TallyObject) :: this + integer(C_INT) :: i + integer(C_INT32_T) :: filt + interface + function tally_get_score_c(tally, i) result(filt) bind(C) + import C_PTR, C_INT + type(C_PTR), value :: tally + integer(C_INT), value :: i + integer(C_INT) :: filt + end function + end interface + filt = tally_get_score_c(this % ptr, i-1) + end function + function tally_get_n_filters(this) result(n) class(TallyObject) :: this integer(C_INT) :: n @@ -721,31 +764,6 @@ contains end function openmc_tally_get_n_realizations - 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 - type(C_PTR), intent(out) :: scores - 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 % score_bins)) then - scores = C_LOC(t % score_bins(1)) - n = size(t % score_bins) - err = 0 - else - err = E_ALLOCATE - call set_errmsg("Tally scores 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_scores - - function openmc_tally_reset(index) result(err) bind(C) ! Reset tally results and number of realizations integer(C_INT32_T), intent(in), value :: index @@ -904,180 +922,6 @@ contains end function openmc_tally_set_nuclides - function openmc_tally_set_scores(index, n, scores) result(err) bind(C) - ! Sets the scores in the tally - integer(C_INT32_T), value :: index - integer(C_INT), value :: n - type(C_PTR), intent(in) :: scores(n) - integer(C_INT) :: err - - integer :: i - integer :: MT - character(C_CHAR), pointer :: string(:) - character(len=:, kind=C_CHAR), allocatable :: score_ - logical :: depletion_rx - - err = E_UNASSIGNED - depletion_rx = .false. - if (index >= 1 .and. index <= size(tallies)) then - associate (t => tallies(index) % obj) - if (allocated(t % score_bins)) deallocate(t % score_bins) - allocate(t % score_bins(n)) - t % n_score_bins = n - - do i = 1, n - ! Convert C string to Fortran string - call c_f_pointer(scores(i), string, [20]) - score_ = to_lower(to_f_string(string)) - - select case (score_) - case ('flux') - t % score_bins(i) = SCORE_FLUX - case ('total', '(n,total)') - t % score_bins(i) = SCORE_TOTAL - case ('scatter') - t % score_bins(i) = SCORE_SCATTER - case ('nu-scatter') - t % score_bins(i) = SCORE_NU_SCATTER - case ('(n,2n)') - t % score_bins(i) = N_2N - depletion_rx = .true. - case ('(n,3n)') - t % score_bins(i) = N_3N - depletion_rx = .true. - case ('(n,4n)') - t % score_bins(i) = N_4N - depletion_rx = .true. - case ('absorption') - t % score_bins(i) = SCORE_ABSORPTION - case ('fission', '18') - t % score_bins(i) = SCORE_FISSION - case ('nu-fission') - t % score_bins(i) = SCORE_NU_FISSION - case ('decay-rate') - t % score_bins(i) = SCORE_DECAY_RATE - case ('delayed-nu-fission') - t % score_bins(i) = SCORE_DELAYED_NU_FISSION - case ('prompt-nu-fission') - t % score_bins(i) = SCORE_PROMPT_NU_FISSION - case ('kappa-fission') - t % score_bins(i) = SCORE_KAPPA_FISSION - case ('inverse-velocity') - t % score_bins(i) = SCORE_INVERSE_VELOCITY - case ('fission-q-prompt') - t % score_bins(i) = SCORE_FISS_Q_PROMPT - case ('fission-q-recoverable') - t % score_bins(i) = SCORE_FISS_Q_RECOV - case ('current') - t % score_bins(i) = SCORE_CURRENT - case ('events') - t % score_bins(i) = SCORE_EVENTS - case ('elastic', '(n,elastic)') - t % score_bins(i) = ELASTIC - case ('(n,2nd)') - t % score_bins(i) = N_2ND - case ('(n,na)') - t % score_bins(i) = N_2NA - case ('(n,n3a)') - t % score_bins(i) = N_N3A - case ('(n,2na)') - t % score_bins(i) = N_2NA - case ('(n,3na)') - t % score_bins(i) = N_3NA - case ('(n,np)') - t % score_bins(i) = N_NP - case ('(n,n2a)') - t % score_bins(i) = N_N2A - case ('(n,2n2a)') - t % score_bins(i) = N_2N2A - case ('(n,nd)') - t % score_bins(i) = N_ND - case ('(n,nt)') - t % score_bins(i) = N_NT - case ('(n,nHe-3)') - t % score_bins(i) = N_N3HE - case ('(n,nd2a)') - t % score_bins(i) = N_ND2A - case ('(n,nt2a)') - t % score_bins(i) = N_NT2A - case ('(n,3nf)') - t % score_bins(i) = N_3NF - case ('(n,2np)') - t % score_bins(i) = N_2NP - case ('(n,3np)') - t % score_bins(i) = N_3NP - case ('(n,n2p)') - t % score_bins(i) = N_N2P - case ('(n,npa)') - t % score_bins(i) = N_NPA - case ('(n,n1)') - t % score_bins(i) = N_N1 - case ('(n,nc)') - t % score_bins(i) = N_NC - case ('(n,gamma)') - t % score_bins(i) = N_GAMMA - depletion_rx = .true. - case ('(n,p)') - t % score_bins(i) = N_P - depletion_rx = .true. - case ('(n,d)') - t % score_bins(i) = N_D - case ('(n,t)') - t % score_bins(i) = N_T - case ('(n,3He)') - t % score_bins(i) = N_3HE - case ('(n,a)') - t % score_bins(i) = N_A - depletion_rx = .true. - case ('(n,2a)') - t % score_bins(i) = N_2A - case ('(n,3a)') - t % score_bins(i) = N_3A - case ('(n,2p)') - t % score_bins(i) = N_2P - case ('(n,pa)') - t % score_bins(i) = N_PA - case ('(n,t2a)') - t % score_bins(i) = N_T2A - case ('(n,d2a)') - t % score_bins(i) = N_D2A - case ('(n,pd)') - t % score_bins(i) = N_PD - case ('(n,pt)') - t % score_bins(i) = N_PT - case ('(n,da)') - t % score_bins(i) = N_DA - case default - ! Assume that user has specified an MT number - MT = int(str_to_int(score_)) - - if (MT /= ERROR_INT) then - ! Specified score was an integer - if (MT > 1) then - t % score_bins(i) = MT - else - err = E_INVALID_ARGUMENT - call set_errmsg("Negative MT number cannot be used as a score.") - end if - - else - err = E_INVALID_ARGUMENT - call set_errmsg("Unknown score: " // trim(score_) // ".") - end if - - end select - end do - - err = 0 - t % depletion_rx = depletion_rx - 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_set_scores - - subroutine openmc_get_tally_next_id(id) bind(C) ! Returns an ID number that has not been used by any other tallies. integer(C_INT32_T), intent(out) :: id diff --git a/src/tallies/trigger.cpp b/src/tallies/trigger.cpp index d8f3ed4b7a..ac69545c61 100644 --- a/src/tallies/trigger.cpp +++ b/src/tallies/trigger.cpp @@ -97,11 +97,8 @@ check_tally_triggers(double& ratio, int& tally_id, int& score) // If this is the most uncertain value, set the output variables. if (this_ratio > ratio) { ratio = this_ratio; - int* scores; - int junk; - err = openmc_tally_get_scores(i_tally, &scores, &junk); - score = scores[trigger.score_index]; - err = openmc_tally_get_id(i_tally, &tally_id); + score = t.scores_[trigger.score_index]; + tally_id = t.id_; } } }