Merge pull request #1050 from shikhar413/cmfd-helper

Added C API bindings and bug fix related to CMFD tallies
This commit is contained in:
Paul Romano 2018-08-20 07:52:55 -05:00 committed by GitHub
commit 6d5ab25a2c
No known key found for this signature in database
GPG key ID: 4AEE18F83AFDEB23
8 changed files with 203 additions and 39 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);
@ -94,15 +95,19 @@ extern "C" {
int openmc_sphharm_filter_set_order(int32_t index, int order);
int openmc_sphharm_filter_set_cosine(int32_t index, const char cosine[]);
int openmc_statepoint_write(const char filename[]);
int openmc_tally_allocate(int32_t index, const char* type);
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);

View file

@ -26,9 +26,15 @@ _dll.openmc_get_tally_index.errcheck = _error_handler
_dll.openmc_global_tallies.argtypes = [POINTER(POINTER(c_double))]
_dll.openmc_global_tallies.restype = c_int
_dll.openmc_global_tallies.errcheck = _error_handler
_dll.openmc_tally_allocate.argtypes = [c_int32, c_char_p]
_dll.openmc_tally_allocate.restype = c_int
_dll.openmc_tally_allocate.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 +53,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 +69,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
@ -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()
@ -170,7 +195,7 @@ class Tally(_FortranObjectWithID):
index = c_int32()
_dll.openmc_extend_tallies(1, index, None)
_dll.openmc_tally_set_type(index, b'generic')
_dll.openmc_tally_allocate(index, b'generic')
index = index.value
else:
index = mapping[uid]._index
@ -190,6 +215,26 @@ 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_set_type(self._index, type.encode())
@property
def estimator(self):
estimator = c_int32()
_dll.openmc_tally_get_estimator(self._index, estimator)
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 +347,10 @@ class Tally(_FortranObjectWithID):
return std_dev
def reset(self):
"""Reset results and num_realizations of tally"""
_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

@ -25,7 +25,7 @@ module openmc_api
use tally_header
use tally_filter_header
use tally_filter
use tally, only: openmc_tally_set_type
use tally, only: openmc_tally_allocate
use simulation
use string, only: to_f_string
use timer_header
@ -60,6 +60,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
@ -83,12 +84,16 @@ module openmc_api
public :: openmc_simulation_init
public :: openmc_source_bank
public :: openmc_source_set_strength
public :: openmc_tally_allocate
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

View file

@ -243,7 +243,7 @@ contains
use error, only: fatal_error, warning
use mesh_header, only: RegularMesh, openmc_extend_meshes
use string
use tally, only: openmc_tally_set_type
use tally, only: openmc_tally_allocate
use tally_header, only: openmc_extend_tallies
use tally_filter_header
use tally_filter
@ -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(:)
@ -433,16 +434,13 @@ contains
! Begin loop around tallies
do i = 1, size(cmfd_tallies)
! Allocate tally
err = openmc_tally_set_type(i_start + i - 1, C_CHAR_'generic' // C_NULL_CHAR)
err = openmc_tally_allocate(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_set_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_set_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_set_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_set_type(i_start + i - 1, C_CHAR_'volume' // C_NULL_CHAR)
! Allocate and set filters
n_filter = 2

View file

@ -1918,6 +1918,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
@ -2116,7 +2117,7 @@ contains
READ_TALLIES: do i = 1, n
! Allocate tally
err = openmc_tally_set_type(i_start + i - 1, &
err = openmc_tally_allocate(i_start + i - 1, &
C_CHAR_'generic' // C_NULL_CHAR)
! Get pointer to tally
@ -2125,19 +2126,15 @@ contains
! Get pointer to tally xml node
node_tal = node_tal_list(i)
! Copy tally id
! Copy and set 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)
err = openmc_tally_set_id(i_start + i - 1, tally_id)
if (err /= 0) call fatal_error(to_f_string(openmc_err_msg))
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
call fatal_error("Two or more tallies use the same unique ID: " &
// to_str(t % id))
end if
! Copy tally name
if (check_for_node(node_tal, "name")) &
call get_node_value(node_tal, "name", t % name)
@ -2837,9 +2834,6 @@ contains
end select
end if
! Add tally to dictionary
call tally_dict % set(t % id, i)
end associate
end do READ_TALLIES

View file

@ -3931,7 +3931,7 @@ contains
! C API FUNCTIONS
!===============================================================================
function openmc_tally_set_type(index, type) result(err) bind(C)
function openmc_tally_allocate(index, type) result(err) bind(C)
! Set the type of the tally
integer(C_INT32_T), value, intent(in) :: index
character(kind=C_CHAR), intent(in) :: type(*)
@ -3964,6 +3964,6 @@ contains
err = E_OUT_OF_BOUNDS
call set_errmsg("Index in tallies array is out of bounds.")
end if
end function openmc_tally_set_type
end function openmc_tally_allocate
end module tally

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_set_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_INVALID_ARGUMENT
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
@ -707,10 +777,17 @@ contains
if (index >= 1 .and. index <= n_tallies) then
if (allocated(tallies(index) % obj)) then
tallies(index) % obj % id = id
call tally_dict % set(id, index)
if (tally_dict % has(id)) then
call set_errmsg("Two or more tallies use the same unique ID: " &
// to_str(id))
err = E_INVALID_ID
else
tallies(index) % obj % id = id
call tally_dict % set(id, index)
if (id > largest_tally_id) largest_tally_id = id
err = 0
err = 0
end if
else
err = E_ALLOCATE
call set_errmsg("Tally type has not been set yet.")
@ -942,4 +1019,42 @@ contains
end if
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
id = largest_tally_id + 1
end subroutine openmc_get_tally_next_id
end module tally_header

View file

@ -159,7 +159,6 @@ def test_tally_mapping(capi_init):
def test_tally(capi_init):
t = openmc.capi.tallies[1]
t.id = 1
assert len(t.filters) == 2
assert isinstance(t.filters[0], openmc.capi.MaterialFilter)
assert isinstance(t.filters[1], openmc.capi.EnergyFilter)
@ -184,7 +183,6 @@ def test_tally(capi_init):
assert t.scores == new_scores
t2 = openmc.capi.tallies[2]
t2.id = 2
assert len(t2.filters) == 2
assert isinstance(t2.filters[0], openmc.capi.ZernikeFilter)
assert isinstance(t2.filters[1], openmc.capi.CellFilter)