diff --git a/include/openmc/tallies/tally.h b/include/openmc/tallies/tally.h index 4b0514fd6..0022fc896 100644 --- a/include/openmc/tallies/tally.h +++ b/include/openmc/tallies/tally.h @@ -30,6 +30,8 @@ public: int32_t n_filter_bins() const {return n_filter_bins_;} + int type_ {TALLY_VOLUME}; //!< volume, surface current + //! Event type that contributes to this tally int estimator_ {ESTIMATOR_TRACKLENGTH}; @@ -64,12 +66,12 @@ extern "C" double total_weight; namespace model { extern std::vector> tallies; + extern std::vector active_tallies; extern std::vector active_analog_tallies; extern std::vector active_tracklength_tallies; - extern std::vector active_meshsurf_tallies; extern std::vector active_collision_tallies; - extern std::vector active_tallies; - extern std::vector active_surface_tallies; + //extern std::vector active_meshsurf_tallies; + //extern std::vector active_surface_tallies; } // Threadprivate variables diff --git a/src/input_xml.F90 b/src/input_xml.F90 index 454956b46..aaab3a460 100644 --- a/src/input_xml.F90 +++ b/src/input_xml.F90 @@ -1394,12 +1394,12 @@ contains &in the same tally as normal surface currents") end if - t % type = TALLY_SURFACE + call t % set_type(TALLY_SURFACE) t % score_bins(j) = SCORE_CURRENT else if (has_meshsurface) then t % score_bins(j) = SCORE_CURRENT - t % type = TALLY_MESH_SURFACE + call t % set_type(TALLY_MESH_SURFACE) ! Check to make sure that current is the only desired response ! for this tally diff --git a/src/state_point.F90 b/src/state_point.F90 index 4fbc09795..85df5b6df 100644 --- a/src/state_point.F90 +++ b/src/state_point.F90 @@ -346,7 +346,7 @@ contains call write_dataset(file_id, "global_tallies", global_tallies) ! Write tallies - if (active_tallies % size() > 0) then + if (active_tallies_size() > 0) then ! Indicate that tallies are on call write_attribute(file_id, "tallies_present", 1) diff --git a/src/tallies/tally.F90 b/src/tallies/tally.F90 index 8420c02cf..226491c0e 100644 --- a/src/tallies/tally.F90 +++ b/src/tallies/tally.F90 @@ -2066,9 +2066,9 @@ contains ! 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_analog_tallies % size() + TALLY_LOOP: do i = 1, active_analog_tallies_size() ! Get index of tally and pointer to tally - i_tally = active_analog_tallies % data(i) + i_tally = active_analog_tallies_data(i) associate (t => tallies(i_tally) % obj) ! Find all valid bins in each filter if they have not already been found @@ -2212,9 +2212,9 @@ contains ! 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_analog_tallies % size() + TALLY_LOOP: do i = 1, active_analog_tallies_size() ! Get index of tally and pointer to tally - i_tally = active_analog_tallies % data(i) + i_tally = active_analog_tallies_data(i) associate (t => tallies(i_tally) % obj) ! Find all valid bins in each filter if they have not already been found @@ -2577,9 +2577,9 @@ contains ! 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_tracklength_tallies % size() + TALLY_LOOP: do i = 1, active_tracklength_tallies_size() ! Get index of tally and pointer to tally - i_tally = active_tracklength_tallies % data(i) + i_tally = active_tracklength_tallies_data(i) associate (t => tallies(i_tally) % obj) ! Find all valid bins in each filter if they have not already been found @@ -2734,9 +2734,9 @@ contains ! 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_collision_tallies % size() + TALLY_LOOP: do i = 1, active_collision_tallies_size() ! Get index of tally and pointer to tally - i_tally = active_collision_tallies % data(i) + i_tally = active_collision_tallies_data(i) associate (t => tallies(i_tally) % obj) ! Find all valid bins in each filter if they have not already been found @@ -3774,8 +3774,8 @@ contains end if ! Accumulate results for each tally - do i = 1, active_tallies % size() - call tallies(active_tallies % data(i)) % obj % accumulate() + do i = 1, active_tallies_size() + call tallies(active_tallies_data(i)) % obj % accumulate() end do end subroutine accumulate_tallies @@ -3790,10 +3790,13 @@ contains integer(C_INT) :: err logical(C_BOOL) :: active - call active_tallies % clear() - call active_analog_tallies % clear() - call active_collision_tallies % clear() - call active_tracklength_tallies % clear() + interface + subroutine setup_active_tallies_c() bind(C) + end subroutine + end interface + + call setup_active_tallies_c() + call active_surface_tallies % clear() call active_meshsurf_tallies % clear() @@ -3801,21 +3804,10 @@ contains associate (t => tallies(i) % obj) err = openmc_tally_get_active(i, active) if (active) then - ! Add tally to active tallies - call active_tallies % push_back(i) - ! Check what type of tally this is and add it to the appropriate list - if (t % type == TALLY_VOLUME) then - if (t % estimator() == ESTIMATOR_ANALOG) then - call active_analog_tallies % push_back(i) - elseif (t % estimator() == ESTIMATOR_TRACKLENGTH) then - call active_tracklength_tallies % push_back(i) - elseif (t % estimator() == ESTIMATOR_COLLISION) then - call active_collision_tallies % push_back(i) - end if - elseif (t % type == TALLY_MESH_SURFACE) then + if (t % type() == TALLY_MESH_SURFACE) then call active_meshsurf_tallies % push_back(i) - elseif (t % type == TALLY_SURFACE) then + elseif (t % type() == TALLY_SURFACE) then call active_surface_tallies % push_back(i) end if diff --git a/src/tallies/tally.cpp b/src/tallies/tally.cpp index cbc7b92f0..a93de78b3 100644 --- a/src/tallies/tally.cpp +++ b/src/tallies/tally.cpp @@ -27,12 +27,12 @@ namespace openmc { namespace model { std::vector> tallies; + std::vector active_tallies; std::vector active_analog_tallies; std::vector active_tracklength_tallies; - std::vector active_meshsurf_tallies; std::vector active_collision_tallies; - std::vector active_tallies; - std::vector active_surface_tallies; + //std::vector active_meshsurf_tallies; + //std::vector active_surface_tallies; } double global_tally_absorption; @@ -177,6 +177,48 @@ void reduce_tally_results() } #endif +extern "C" void +setup_active_tallies_c() +{ + model::active_tallies.clear(); + model::active_analog_tallies.clear(); + model::active_tracklength_tallies.clear(); + model::active_collision_tallies.clear(); + //model::active_meshsurf_tallies.clear(); + //model::active_surface_tallies.clear(); + + //TODO: off-by-one all through here + for (auto i = 0; i < model::tallies.size(); ++i) { + const auto& tally {*model::tallies[i]}; + + if (tally.active_) { + model::active_tallies.push_back(i + 1); + switch (tally.type_) { + + case TALLY_VOLUME: + switch (tally.estimator_) { + case ESTIMATOR_ANALOG: + model::active_analog_tallies.push_back(i + 1); + break; + case ESTIMATOR_TRACKLENGTH: + model::active_tracklength_tallies.push_back(i + 1); + break; + case ESTIMATOR_COLLISION: + model::active_collision_tallies.push_back(i + 1); + } + break; + + //case TALLY_MESH_SURFACE: + // model::active_meshsurf_tallies.push_back(i + 1); + // break; + + //case TALLY_SURFACE: + // model::active_surface_tallies.push_back(i + 1); + } + } + } +} + extern "C" void free_memory_tally_c() { @@ -189,12 +231,51 @@ free_memory_tally_c() model::tally_filters.clear(); model::tallies.clear(); + + model::active_tallies.clear(); + model::active_analog_tallies.clear(); + model::active_tracklength_tallies.clear(); + model::active_collision_tallies.clear(); + //model::active_meshsurf_tallies.clear(); + //model::active_surface_tallies.clear(); } //============================================================================== // C-API functions //============================================================================== +extern "C" int +openmc_tally_get_type(int32_t index, int32_t* type) +{ + 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 + *type = model::tallies[index-1]->type_; +} + +extern "C" int +openmc_tally_set_type(int32_t index, const char* type) +{ + if (index < 1 || index > model::tallies.size()) { + set_errmsg("Index in tallies array is out of bounds."); + return OPENMC_E_OUT_OF_BOUNDS; + } + if (strcmp(type, "volume") == 0) { + model::tallies[index-1]->type_ = TALLY_VOLUME; + } else if (strcmp(type, "mesh-surface") == 0) { + model::tallies[index-1]->type_ = TALLY_MESH_SURFACE; + } else if (strcmp(type, "surface") == 0) { + model::tallies[index-1]->type_ = TALLY_SURFACE; + } else { + std::stringstream errmsg; + errmsg << "Unknown tally type: " << type; + set_errmsg(errmsg); + return OPENMC_E_INVALID_ARGUMENT; + } +} + extern "C" int openmc_tally_get_active(int32_t index, bool* active) { @@ -266,6 +347,34 @@ extern "C" { model::tallies.push_back(std::make_unique()); } + int active_tallies_data(int i) + {return model::active_tallies[i-1];} + + int active_tallies_size() + {return model::active_tallies.size();} + + int active_analog_tallies_data(int i) + {return model::active_analog_tallies[i-1];} + + int active_analog_tallies_size() + {return model::active_analog_tallies.size();} + + int active_tracklength_tallies_data(int i) + {return model::active_tracklength_tallies[i-1];} + + int active_tracklength_tallies_size() + {return model::active_tracklength_tallies.size();} + + int active_collision_tallies_data(int i) + {return model::active_collision_tallies[i-1];} + + int active_collision_tallies_size() + {return model::active_collision_tallies.size();} + + int tally_get_type_c(Tally* tally) {return tally->type_;} + + void tally_set_type_c(Tally* tally, int type) {tally->type_ = type;} + int tally_get_estimator_c(Tally* tally) {return tally->estimator_;} void tally_set_estimator_c(Tally* tally, int e) {tally->estimator_ = e;} diff --git a/src/tallies/tally_header.F90 b/src/tallies/tally_header.F90 index aa929e8c4..2e57d243b 100644 --- a/src/tallies/tally_header.F90 +++ b/src/tallies/tally_header.F90 @@ -16,29 +16,6 @@ module tally_header use trigger_header, only: TriggerObject implicit none - private - public :: allocate_tally_results - public :: free_memory_tally - public :: openmc_extend_tallies - public :: openmc_get_tally_index - public :: openmc_get_tally_next_id - public :: openmc_global_tallies - public :: openmc_tally_get_active - public :: openmc_tally_get_estimator - public :: openmc_tally_get_id - public :: openmc_tally_get_n_realizations - public :: openmc_tally_get_nuclides - public :: openmc_tally_get_scores - public :: openmc_tally_get_type - public :: openmc_tally_reset - public :: openmc_tally_results - public :: openmc_tally_set_active - public :: openmc_tally_set_estimator - public :: openmc_tally_set_filters - public :: openmc_tally_set_id - public :: openmc_tally_set_nuclides - public :: openmc_tally_set_scores - public :: openmc_tally_set_type interface function openmc_tally_set_filters(index, n, filter_indices) result(err) bind(C) @@ -62,6 +39,50 @@ module tally_header logical(C_BOOL), intent(out) :: active integer(C_INT) :: err end function + + function active_tallies_size() result(size) bind(C) + import C_INT + integer(C_INT) :: size + end function + + function active_tallies_data(i) result(tally) bind(C) + import C_INT + integer(C_INT), value :: i + integer(C_INT) :: tally + end function + + function active_analog_tallies_size() result(size) bind(C) + import C_INT + integer(C_INT) :: size + end function + + function active_analog_tallies_data(i) result(tally) bind(C) + import C_INT + integer(C_INT), value :: i + integer(C_INT) :: tally + end function + + function active_tracklength_tallies_size() result(size) bind(C) + import C_INT + integer(C_INT) :: size + end function + + function active_tracklength_tallies_data(i) result(tally) bind(C) + import C_INT + integer(C_INT), value :: i + integer(C_INT) :: tally + end function + + function active_collision_tallies_size() result(size) bind(C) + import C_INT + integer(C_INT) :: size + end function + + function active_collision_tallies_data(i) result(tally) bind(C) + import C_INT + integer(C_INT), value :: i + integer(C_INT) :: tally + end function end interface !=============================================================================== @@ -77,10 +98,7 @@ module tally_header integer :: id ! user-defined identifier character(len=104) :: name = "" ! user-defined name - integer :: type = TALLY_VOLUME ! volume, surface current - !integer :: estimator = ESTIMATOR_TRACKLENGTH ! collision, track-length real(8) :: volume ! volume of region - !logical :: active = .false. logical :: depletion_rx = .false. ! has depletion reactions, e.g. (n,2n) ! Individual nuclides to tally @@ -112,6 +130,8 @@ module tally_header procedure :: allocate_results => tally_allocate_results procedure :: read_results_hdf5 => tally_read_results_hdf5 procedure :: write_results_hdf5 => tally_write_results_hdf5 + procedure :: type => tally_get_type + procedure :: set_type => tally_set_type procedure :: estimator => tally_get_estimator procedure :: set_estimator => tally_set_estimator procedure :: n_filters => tally_get_n_filters @@ -163,11 +183,7 @@ module tally_header !$omp& global_tally_tracklength, global_tally_leakage) ! Active tally lists - type(VectorInt), public :: active_analog_tallies - type(VectorInt), public :: active_tracklength_tallies type(VectorInt), public :: active_meshsurf_tallies - type(VectorInt), public :: active_collision_tallies - type(VectorInt), public :: active_tallies type(VectorInt), public :: active_surface_tallies ! Normalization for statistics @@ -298,6 +314,32 @@ contains end subroutine tally_allocate_results + function tally_get_type(this) result(t) + class(TallyObject) :: this + integer(C_INT) :: t + interface + function tally_get_type_c(tally) result(t) bind(C) + import C_PTR, C_INT + type(C_PTR), value :: tally + integer(C_INT) :: t + end function + end interface + t = tally_get_type_c(this % ptr) + end function + + subroutine tally_set_type(this, t) + class(TallyObject) :: this + integer(C_INT) :: t + interface + subroutine tally_set_type_c(tally, t) bind(C) + import C_PTR, C_INT + type(C_PTR), value :: tally + integer(C_INT), value :: t + end subroutine + end interface + call tally_set_type_c(this % ptr, t) + end subroutine + function tally_get_estimator(this) result(e) class(TallyObject) :: this integer(C_INT) :: e @@ -513,12 +555,8 @@ contains if (allocated(global_tallies)) deallocate(global_tallies) ! Deallocate tally node lists - call active_analog_tallies % clear() - call active_tracklength_tallies % clear() call active_meshsurf_tallies % clear() - call active_collision_tallies % clear() call active_surface_tallies % clear() - call active_tallies % clear() end subroutine free_memory_tally !=============================================================================== @@ -701,22 +739,6 @@ contains end function openmc_tally_get_scores - function openmc_tally_get_type(index, type) result(err) bind(C) - ! Return the type of a tally - integer(C_INT32_T), value :: index - integer(C_INT32_T), intent(out) :: type - integer(C_INT) :: err - - if (index >= 1 .and. index <= size(tallies)) then - type = tallies(index) % obj % type - err = 0 - else - err = E_OUT_OF_BOUNDS - call set_errmsg('Index in tallies array is out of bounds.') - end if - end function openmc_tally_get_type - - function openmc_tally_reset(index) result(err) bind(C) ! Reset tally results and number of realizations integer(C_INT32_T), intent(in), value :: index @@ -1048,37 +1070,6 @@ contains end function openmc_tally_set_scores - function openmc_tally_set_type(index, type) result(err) bind(C) - ! Update the type of a tally that is already allocated - integer(C_INT32_T), value, intent(in) :: index - character(kind=C_CHAR), intent(in) :: type(*) - integer(C_INT) :: err - - character(:), allocatable :: type_ - - ! Convert C string to Fortran string - type_ = to_f_string(type) - - err = 0 - if (index >= 1 .and. index <= size(tallies)) then - select case (type_) - case ('volume') - tallies(index) % obj % type = TALLY_VOLUME - case ('mesh-surface') - tallies(index) % obj % type = TALLY_MESH_SURFACE - case ('surface') - tallies(index) % obj % type = TALLY_SURFACE - case default - err = E_INVALID_ARGUMENT - call set_errmsg("Unknown tally type: " // trim(type_)) - end select - else - err = E_OUT_OF_BOUNDS - call set_errmsg("Index in tally array is out of bounds.") - end if - end function openmc_tally_set_type - - 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.F90 b/src/tallies/trigger.F90 index 31281a64f..4c1efcc2e 100644 --- a/src/tallies/trigger.F90 +++ b/src/tallies/trigger.F90 @@ -166,7 +166,7 @@ contains trigger % variance = ZERO ! Mesh current tally triggers require special treatment - if (t % type == TALLY_MESH_SURFACE) then + if (t % type() == TALLY_MESH_SURFACE) then call compute_tally_current(t, trigger) else diff --git a/src/tracking.F90 b/src/tracking.F90 index bf8332852..f540687f3 100644 --- a/src/tracking.F90 +++ b/src/tracking.F90 @@ -86,7 +86,7 @@ contains endif ! Every particle starts with no accumulated flux derivative. - if (active_tallies % size() > 0) call zero_flux_derivs() + if (active_tallies_size() > 0) call zero_flux_derivs() EVENT_LOOP: do ! Set the random number stream @@ -172,7 +172,7 @@ contains end do ! Score track-length tallies - if (active_tracklength_tallies % size() > 0) then + if (active_tracklength_tallies_size() > 0) then call score_tracklength_tally(p, distance) end if @@ -183,7 +183,7 @@ contains end if ! Score flux derivative accumulators for differential tallies. - if (active_tallies % size() > 0) call score_track_derivative(p, distance) + if (active_tallies_size() > 0) call score_track_derivative(p, distance) if (d_collision > d_boundary) then ! ==================================================================== @@ -241,8 +241,8 @@ contains ! Score collision estimator tallies -- this is done after a collision ! has occurred rather than before because we need information on the ! outgoing energy for any tallies with an outgoing energy filter - if (active_collision_tallies % size() > 0) call score_collision_tally(p) - if (active_analog_tallies % size() > 0) call score_analog_tally(p) + if (active_collision_tallies_size() > 0) call score_collision_tally(p) + if (active_analog_tallies_size() > 0) call score_analog_tally(p) ! Reset banked weight during collision p % n_bank = 0 @@ -273,7 +273,7 @@ contains end do ! Score flux derivative accumulators for differential tallies. - if (active_tallies % size() > 0) call score_collision_derivative(p) + if (active_tallies_size() > 0) call score_collision_derivative(p) end if ! If particle has too many events, display warning and kill it diff --git a/tests/unit_tests/test_capi.py b/tests/unit_tests/test_capi.py index f5185e23f..42c0639cb 100644 --- a/tests/unit_tests/test_capi.py +++ b/tests/unit_tests/test_capi.py @@ -162,6 +162,7 @@ def test_tally_mapping(capi_init): def test_tally(capi_init): t = openmc.capi.tallies[1] + assert t.type == 'volume' assert len(t.filters) == 2 assert isinstance(t.filters[0], openmc.capi.MaterialFilter) assert isinstance(t.filters[1], openmc.capi.EnergyFilter)