Extract code from shikhar413::cmfd-capi relevant to PR

This commit is contained in:
Shikhar Kumar 2018-08-16 21:45:47 -04:00
parent c56b91bc77
commit 1ec0a806a0
6 changed files with 205 additions and 34 deletions

View file

@ -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);

View file

@ -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

View file

@ -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

View file

@ -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

View file

@ -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 <scores> 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

View file

@ -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