From 1ec0a806a0b0bba223d64735458fab9fe6a48248 Mon Sep 17 00:00:00 2001 From: Shikhar Kumar Date: Thu, 16 Aug 2018 21:45:47 -0400 Subject: [PATCH] Extract code from shikhar413::cmfd-capi relevant to PR --- include/openmc.h | 5 ++ openmc/capi/tally.py | 52 ++++++++++++++++ src/api.F90 | 5 ++ src/cmfd_input.F90 | 24 ++++---- src/input_xml.F90 | 37 +++++------ src/tallies/tally_header.F90 | 116 ++++++++++++++++++++++++++++++++++- 6 files changed, 205 insertions(+), 34 deletions(-) diff --git a/include/openmc.h b/include/openmc.h index 1f919d7c61..36bfd81cbf 100644 --- a/include/openmc.h +++ b/include/openmc.h @@ -47,6 +47,7 @@ extern "C" { int openmc_get_nuclide_index(const char name[], int* index); int64_t openmc_get_seed(); int openmc_get_tally_index(int32_t id, int32_t* index); + void openmc_get_tally_next_id(int32_t* id); int openmc_hard_reset(); int openmc_init(int argc, char* argv[], const void* intracomm); int openmc_init_f(const int* intracomm); @@ -93,19 +94,23 @@ extern "C" { int openmc_sphharm_filter_set_cosine(int32_t index, const char cosine[]); int openmc_statepoint_write(const char filename[]); int openmc_tally_get_active(int32_t index, bool* active); + int openmc_tally_get_estimator(int32_t index, int32_t* estimator); int openmc_tally_get_id(int32_t index, int32_t* id); int openmc_tally_get_filters(int32_t index, int32_t** indices, int* n); int openmc_tally_get_n_realizations(int32_t index, int32_t* n); int openmc_tally_get_nuclides(int32_t index, int** nuclides, int* n); int openmc_tally_get_scores(int32_t index, int** scores, int* n); + int openmc_tally_get_type(int32_t index, int32_t* type); int openmc_tally_reset(int32_t index); int openmc_tally_results(int32_t index, double** ptr, int shape_[3]); int openmc_tally_set_active(int32_t index, bool active); + int openmc_tally_set_estimator(int32_t index, const char* estimator); int openmc_tally_set_filters(int32_t index, int n, const int32_t* indices); int openmc_tally_set_id(int32_t index, int32_t id); int openmc_tally_set_nuclides(int32_t index, int n, const char** nuclides); int openmc_tally_set_scores(int32_t index, int n, const char** scores); int openmc_tally_set_type(int32_t index, const char* type); + int openmc_tally_update_type(int32_t index, const char* type); int openmc_zernike_filter_get_order(int32_t index, int* order); int openmc_zernike_filter_get_params(int32_t index, double* x, double* y, double* r); int openmc_zernike_filter_set_order(int32_t index, int order); diff --git a/openmc/capi/tally.py b/openmc/capi/tally.py index d4d70af718..e494e33166 100644 --- a/openmc/capi/tally.py +++ b/openmc/capi/tally.py @@ -29,6 +29,9 @@ _dll.openmc_global_tallies.errcheck = _error_handler _dll.openmc_tally_get_active.argtypes = [c_int32, POINTER(c_bool)] _dll.openmc_tally_get_active.restype = c_int _dll.openmc_tally_get_active.errcheck = _error_handler +_dll.openmc_tally_get_estimator.argtypes = [c_int32, POINTER(c_int32)] +_dll.openmc_tally_get_estimator.restype = c_int +_dll.openmc_tally_get_estimator.errcheck = _error_handler _dll.openmc_tally_get_id.argtypes = [c_int32, POINTER(c_int32)] _dll.openmc_tally_get_id.restype = c_int _dll.openmc_tally_get_id.errcheck = _error_handler @@ -47,6 +50,12 @@ _dll.openmc_tally_get_scores.argtypes = [ c_int32, POINTER(POINTER(c_int)), POINTER(c_int)] _dll.openmc_tally_get_scores.restype = c_int _dll.openmc_tally_get_scores.errcheck = _error_handler +_dll.openmc_tally_get_type.argtypes = [c_int32, POINTER(c_int32)] +_dll.openmc_tally_get_type.restype = c_int +_dll.openmc_tally_get_type.errcheck = _error_handler +_dll.openmc_tally_reset.argtypes = [c_int32] +_dll.openmc_tally_reset.restype = c_int +_dll.openmc_tally_reset.errcheck = _error_handler _dll.openmc_tally_results.argtypes = [ c_int32, POINTER(POINTER(c_double)), POINTER(c_int*3)] _dll.openmc_tally_results.restype = c_int @@ -57,6 +66,9 @@ _dll.openmc_tally_set_active.errcheck = _error_handler _dll.openmc_tally_set_filters.argtypes = [c_int32, c_int, POINTER(c_int32)] _dll.openmc_tally_set_filters.restype = c_int _dll.openmc_tally_set_filters.errcheck = _error_handler +_dll.openmc_tally_set_estimator.argtypes = [c_int32, c_char_p] +_dll.openmc_tally_set_estimator.restype = c_int +_dll.openmc_tally_set_estimator.errcheck = _error_handler _dll.openmc_tally_set_id.argtypes = [c_int32, c_int32] _dll.openmc_tally_set_id.restype = c_int _dll.openmc_tally_set_id.errcheck = _error_handler @@ -69,6 +81,9 @@ _dll.openmc_tally_set_scores.errcheck = _error_handler _dll.openmc_tally_set_type.argtypes = [c_int32, c_char_p] _dll.openmc_tally_set_type.restype = c_int _dll.openmc_tally_set_type.errcheck = _error_handler +_dll.openmc_tally_update_type.argtypes = [c_int32, c_char_p] +_dll.openmc_tally_update_type.restype = c_int +_dll.openmc_tally_update_type.errcheck = _error_handler _SCORES = { @@ -78,6 +93,12 @@ _SCORES = { -12: 'prompt-nu-fission', -13: 'inverse-velocity', -14: 'fission-q-prompt', -15: 'fission-q-recoverable', -16: 'decay-rate' } +_ESTIMATORS = { + 1: 'analog', 2: 'tracklength', 3: 'collision' +} +_TALLY_TYPES = { + 1: 'volume', 2: 'mesh-surface', 3: 'surface' +} def global_tallies(): @@ -140,6 +161,8 @@ class Tally(_FortranObjectWithID): ---------- id : int ID of the tally + estimator: str + Estimator type of tally (analog, tracklength, collision) filters : list List of tally filters mean : numpy.ndarray @@ -152,6 +175,8 @@ class Tally(_FortranObjectWithID): Array of tally results std_dev : numpy.ndarray An array containing the sample standard deviation for each bin + type : str + Type of tally (volume, mesh_surface, surface) """ __instances = WeakValueDictionary() @@ -190,6 +215,30 @@ class Tally(_FortranObjectWithID): _dll.openmc_tally_get_active(self._index, active) return active.value + @property + def type(self): + type = c_int32() + _dll.openmc_tally_get_type(self._index, type) + return _TALLY_TYPES[type.value] + + @type.setter + def type(self, type): + _dll.openmc_tally_update_type(self._index, type.encode()) + + @property + def estimator(self): + estimator = c_int32() + try: + _dll.openmc_tally_get_estimator(self._index, estimator) + except AllocationError: + return "" + else: + return _ESTIMATORS[estimator.value] + + @estimator.setter + def estimator(self, estimator): + _dll.openmc_tally_set_estimator(self._index, estimator.encode()) + @active.setter def active(self, active): _dll.openmc_tally_set_active(self._index, active) @@ -302,6 +351,9 @@ class Tally(_FortranObjectWithID): return std_dev + def reset(self): + _dll.openmc_tally_reset(self._index) + def ci_width(self, alpha=0.05): """Confidence interval half-width based on a Student t distribution diff --git a/src/api.F90 b/src/api.F90 index f7e896b7c8..61cc6150c3 100644 --- a/src/api.F90 +++ b/src/api.F90 @@ -62,6 +62,7 @@ module openmc_api public :: openmc_get_nuclide_index public :: openmc_get_seed public :: openmc_get_tally_index + public :: openmc_get_tally_next_id public :: openmc_global_tallies public :: openmc_hard_reset public :: openmc_init_f @@ -85,17 +86,21 @@ module openmc_api public :: openmc_simulation_init public :: openmc_source_bank public :: openmc_source_set_strength + public :: openmc_tally_get_estimator public :: openmc_tally_get_id public :: openmc_tally_get_filters public :: openmc_tally_get_n_realizations public :: openmc_tally_get_nuclides public :: openmc_tally_get_scores + public :: openmc_tally_get_type public :: openmc_tally_results + 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 + public :: openmc_tally_update_type contains diff --git a/src/cmfd_input.F90 b/src/cmfd_input.F90 index c1965b2206..33cd01a9e9 100644 --- a/src/cmfd_input.F90 +++ b/src/cmfd_input.F90 @@ -262,6 +262,7 @@ contains integer(C_INT) :: err integer :: i_filt ! index in filters array integer :: filt_id + integer :: tally_id integer :: iarray3(3) ! temp integer array real(8) :: rarray3(3) ! temp double array real(C_DOUBLE), allocatable :: energies(:) @@ -434,15 +435,12 @@ contains do i = 1, size(cmfd_tallies) ! Allocate tally err = openmc_tally_set_type(i_start + i - 1, C_CHAR_'generic' // C_NULL_CHAR) + call openmc_get_tally_next_id(tally_id) + err = openmc_tally_set_id(i_start + i - 1, tally_id) ! Point t to tally variable associate (t => cmfd_tallies(i) % obj) - ! Set reset property - if (check_for_node(root, "reset")) then - call get_node_value(root, "reset", t % reset) - end if - ! Set the incoming energy mesh filter index in the tally find_filter ! array n_filter = 1 @@ -464,10 +462,10 @@ contains t % name = "CMFD flux, total" ! Set tally estimator to analog - t % estimator = ESTIMATOR_ANALOG + err = openmc_tally_set_estimator(i_start + i - 1, C_CHAR_'analog' // C_NULL_CHAR) ! Set tally type to volume - t % type = TALLY_VOLUME + err = openmc_tally_update_type(i_start + i - 1, C_CHAR_'volume' // C_NULL_CHAR) ! Allocate and set filters allocate(filter_indices(n_filter)) @@ -492,10 +490,10 @@ contains t % name = "CMFD neutron production" ! Set tally estimator to analog - t % estimator = ESTIMATOR_ANALOG + err = openmc_tally_set_estimator(i_start + i - 1, C_CHAR_'analog' // C_NULL_CHAR) ! Set tally type to volume - t % type = TALLY_VOLUME + err = openmc_tally_update_type(i_start + i - 1, C_CHAR_'volume' // C_NULL_CHAR) ! Set the incoming energy mesh filter index in the tally find_filter ! array @@ -527,7 +525,7 @@ contains t % name = "CMFD surface currents" ! Set tally estimator to analog - t % estimator = ESTIMATOR_ANALOG + err = openmc_tally_set_estimator(i_start + i - 1, C_CHAR_'analog' // C_NULL_CHAR) ! Allocate and set filters allocate(filter_indices(n_filter)) @@ -544,17 +542,17 @@ contains ! Set macro bins t % score_bins(1) = SCORE_CURRENT - t % type = TALLY_MESH_SURFACE + err = openmc_tally_update_type(i_start + i - 1, C_CHAR_'mesh-surface' // C_NULL_CHAR) else if (i == 4) then ! Set name t % name = "CMFD P1 scatter" ! Set tally estimator to analog - t % estimator = ESTIMATOR_ANALOG + err = openmc_tally_set_estimator(i_start + i - 1, C_CHAR_'analog' // C_NULL_CHAR) ! Set tally type to volume - t % type = TALLY_VOLUME + err = openmc_tally_update_type(i_start + i - 1, C_CHAR_'volume' // C_NULL_CHAR) ! Allocate and set filters n_filter = 2 diff --git a/src/input_xml.F90 b/src/input_xml.F90 index da33a6627f..a6db11c422 100644 --- a/src/input_xml.F90 +++ b/src/input_xml.F90 @@ -1960,6 +1960,7 @@ contains integer :: k ! another loop index integer :: l ! loop over bins integer :: filter_id ! user-specified identifier for filter + integer :: tally_id ! user-specified identifier for filter integer :: i_filt ! index in filters array integer :: i_elem ! index of entry in dictionary integer :: n ! size of arrays in mesh specification @@ -2169,15 +2170,15 @@ contains ! Copy tally id if (check_for_node(node_tal, "id")) then - call get_node_value(node_tal, "id", t % id) + call get_node_value(node_tal, "id", tally_id) else call fatal_error("Must specify id for tally in tally XML file.") end if ! Check to make sure 'id' hasn't been used - if (tally_dict % has(t % id)) then + if (tally_dict % has(tally_id)) then call fatal_error("Two or more tallies use the same unique ID: " & - // to_str(t % id)) + // to_str(tally_id)) end if ! Copy tally name @@ -2216,7 +2217,7 @@ contains else call fatal_error("Could not find filter " & // trim(to_str(temp_filter(j))) // " specified on tally " & - // trim(to_str(t % id))) + // trim(to_str(tally_id))) end if ! Store the index of the filter @@ -2270,7 +2271,7 @@ contains if (.not. nuclide_dict % has(word)) then call fatal_error("Could not find the nuclide " & // trim(word) // " specified in tally " & - // trim(to_str(t % id)) // " in any material.") + // trim(to_str(tally_id)) // " in any material.") end if ! Set bin to index in nuclides array @@ -2580,7 +2581,7 @@ contains 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))) + // trim(to_str(tally_id))) end if end do end do @@ -2627,7 +2628,7 @@ contains end if else call fatal_error("No specified on tally " & - // trim(to_str(t % id)) // ".") + // trim(to_str(tally_id)) // ".") end if ! Check for a tally derivative. @@ -2650,7 +2651,7 @@ contains if (j == size(tally_derivs)) then call fatal_error("Could not find derivative " & // trim(to_str(t % deriv)) // " specified on tally " & - // trim(to_str(t % id))) + // trim(to_str(tally_id))) end if end do @@ -2658,7 +2659,7 @@ contains .or. tally_derivs(t % deriv) % variable == DIFF_TEMPERATURE) then if (any(t % nuclide_bins == -1)) then if (t % find_filter(FILTER_ENERGYOUT) > 0) then - call fatal_error("Error on tally " // trim(to_str(t % id)) & + call fatal_error("Error on tally " // trim(to_str(tally_id)) & // ": Cannot use a 'nuclide_density' or 'temperature' & &derivative on a tally with an outgoing energy filter and & &'total' nuclide rate. Instead, tally each nuclide in the & @@ -2735,7 +2736,7 @@ contains temp_str = to_lower(temp_str) else call fatal_error("Must specify trigger type for tally " // & - trim(to_str(t % id)) // " in tally XML file.") + trim(to_str(tally_id)) // " in tally XML file.") end if ! Get the convergence threshold for the trigger @@ -2743,7 +2744,7 @@ contains call get_node_value(node_trigger, "threshold", threshold) else call fatal_error("Must specify trigger threshold for tally " // & - trim(to_str(t % id)) // " in tally XML file.") + trim(to_str(tally_id)) // " in tally XML file.") end if ! Get list scores for this trigger @@ -2787,7 +2788,7 @@ contains t % triggers(trig_ind) % type = RELATIVE_ERROR case default call fatal_error("Unknown trigger type " // & - trim(temp_str) // " in tally " // trim(to_str(t % id))) + trim(temp_str) // " in tally " // trim(to_str(tally_id))) end select ! Store the trigger convergence threshold @@ -2808,7 +2809,7 @@ contains ! Check if an invalid score was set for the trigger if (t % triggers(trig_ind) % score_index == 0) then call fatal_error("The trigger score " // trim(score_name) // & - " is not set for tally " // trim(to_str(t % id))) + " is not set for tally " // trim(to_str(tally_id))) end if ! Store the trigger convergence threshold @@ -2824,7 +2825,7 @@ contains t % triggers(trig_ind) % type = RELATIVE_ERROR case default call fatal_error("Unknown trigger type " // trim(temp_str) // & - " in tally " // trim(to_str(t % id))) + " in tally " // trim(to_str(tally_id))) end select ! Increment the overall trigger index @@ -2856,7 +2857,7 @@ contains ! tally needs post-collision information if (t % estimator == ESTIMATOR_ANALOG) then call fatal_error("Cannot use track-length estimator for tally " & - // to_str(t % id)) + // to_str(tally_id)) end if ! Set estimator to track-length estimator @@ -2867,7 +2868,7 @@ contains ! tally needs post-collision information if (t % estimator == ESTIMATOR_ANALOG) then call fatal_error("Cannot use collision estimator for tally " & - // to_str(t % id)) + // to_str(tally_id)) end if ! Set estimator to collision estimator @@ -2879,8 +2880,8 @@ contains end select end if - ! Add tally to dictionary - call tally_dict % set(t % id, i) + ! Set tally id + err = openmc_tally_set_id(i_start + i - 1, tally_id) end associate end do READ_TALLIES diff --git a/src/tallies/tally_header.F90 b/src/tallies/tally_header.F90 index c64fc9d5ab..d22a53c6a2 100644 --- a/src/tallies/tally_header.F90 +++ b/src/tallies/tally_header.F90 @@ -22,20 +22,25 @@ module tally_header 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_filters 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_update_type !=============================================================================== ! TALLYOBJECT describes a user-specified tally. The region of phase space to @@ -86,9 +91,6 @@ module tally_header integer :: total_score_bins real(C_DOUBLE), allocatable :: results(:,:,:) - ! reset property - allows a tally to be reset after every batch - logical :: reset = .false. - ! Number of realizations of tally random variables integer :: n_realizations = 0 @@ -118,6 +120,10 @@ module tally_header ! Dictionary that maps user IDs to indices in 'tallies' type(DictIntInt), public :: tally_dict + ! The largest tally ID that has been specified in the system. This is useful + ! in case the code needs to find an ID for a new tally. + integer :: largest_tally_id + ! Global tallies ! 1) collision estimate of k-eff ! 2) absorption estimate of k-eff @@ -403,6 +409,7 @@ contains n_tallies = 0 if (allocated(tallies)) deallocate(tallies) call tally_dict % clear() + largest_tally_id = 0 if (allocated(global_tallies)) deallocate(global_tallies) @@ -505,6 +512,22 @@ contains end function openmc_tally_get_active + function openmc_tally_get_estimator(index, estimator) result(err) bind(C) + ! Return the type of estimator of a tally + integer(C_INT32_T), value :: index + integer(C_INT32_T), intent(out) :: estimator + integer(C_INT) :: err + + if (index >= 1 .and. index <= size(tallies)) then + estimator = tallies(index) % obj % estimator + 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_estimator + + function openmc_tally_get_id(index, id) result(err) bind(C) ! Return the ID of a tally integer(C_INT32_T), value :: index @@ -612,6 +635,22 @@ 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 @@ -656,6 +695,37 @@ contains end function openmc_tally_results + function openmc_tally_set_estimator(index, estimator) result(err) bind(C) + ! Set the type of estimator a tally + integer(C_INT32_T), value, intent(in) :: index + character(kind=C_CHAR), intent(in) :: estimator(*) + integer(C_INT) :: err + + character(:), allocatable :: estimator_ + + ! Convert C string to Fortran string + estimator_ = to_f_string(estimator) + + err = 0 + if (index >= 1 .and. index <= size(tallies)) then + select case (estimator_) + case ('analog') + tallies(index) % obj % estimator = ESTIMATOR_ANALOG + case ('tracklength') + tallies(index) % obj % estimator = ESTIMATOR_TRACKLENGTH + case ('collision') + tallies(index) % obj % estimator = ESTIMATOR_COLLISION + case default + err = E_UNASSIGNED + call set_errmsg("Unknown tally estimator: " // trim(estimator_)) + 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_estimator + + function openmc_tally_set_filters(index, n, filter_indices) result(err) bind(C) ! Set the list of filters for a tally integer(C_INT32_T), value, intent(in) :: index @@ -709,6 +779,7 @@ contains if (allocated(tallies(index) % obj)) then tallies(index) % obj % id = id call tally_dict % set(id, index) + if (id > largest_tally_id) largest_tally_id = id err = 0 else @@ -942,4 +1013,43 @@ contains end if end function openmc_tally_set_scores + + function openmc_tally_update_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_UNASSIGNED + 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_update_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 + + id = largest_tally_id + 1 + end subroutine openmc_get_tally_next_id + end module tally_header