crp/OpenMC
1
0
Fork 0
mirror of https://github.com/openmc-dev/openmc.git synced 2026-07-28 14:15:42 -04:00
Code Issues Projects Releases Packages Wiki Activity Actions 5

Make k-effective accessible from C API

Browse source
This commit is contained in:
Paul Romano 2017-07-25 22:13:41 -05:00
parent 35657d5e1b
commit 7bc61c1f41
8 changed files with 54 additions and 26 deletions
Download patch file Download diff file Expand all files Collapse all files

16
openmc/capi.py
View file

@ -43,6 +43,8 @@ class OpenMCLibrary(object):
self._dll.openmc_find.restype = None
self._dll.openmc_init.argtypes = [POINTER(c_int)]
self._dll.openmc_init.restype = None
self._dll.openmc_get_keff.argtypes = [POINTER(c_double*2)]
self._dll.openmc_get_keff.restype = c_int
self._dll.openmc_load_nuclide.argtypes = [c_char_p]
self._dll.openmc_load_nuclide.restype = c_int
self._dll.openmc_material_add_nuclide.argtypes = [
@ -141,7 +143,21 @@ class OpenMCLibrary(object):
else:
return self._dll.openmc_init(None)
def keff(self):
"""Return the calculated k-eigenvalue and its standard deviation.
Returns
-------
tuple
Mean k-eigenvalue and standard deviation of the mean
"""
k = (c_double*2)()
err = self._dll.openmc_get_keff(k)
return tuple(k)
def load_nuclide(self, name):
"""Load cross section data for a nuclide.
Parameters

3
src/api.F90
View file

@ -5,7 +5,7 @@ module openmc_api
use hdf5, only: HID_T
use constants, only: K_BOLTZMANN
use eigenvalue, only: k_sum
use eigenvalue, only: k_sum, openmc_get_keff
use finalize, only: openmc_finalize
use geometry, only: find_cell
use global
@ -23,6 +23,7 @@ module openmc_api
public :: openmc_cell_set_temperature
public :: openmc_finalize
public :: openmc_find
public :: openmc_get_keff
public :: openmc_init
public :: openmc_load_nuclide
public :: openmc_material_add_nuclide

20
src/eigenvalue.F90
View file

@ -1,5 +1,6 @@
module eigenvalue
use, intrinsic :: ISO_C_BINDING
use algorithm, only: binary_search
use constants, only: ZERO
@ -439,8 +440,8 @@ contains
end subroutine calculate_average_keff
!===============================================================================
! CALCULATE_COMBINED_KEFF calculates a minimum variance estimate of k-effective
! based on a linear combination of the collision, absorption, and tracklength
! OPENMC_GET_KEFF calculates a minimum variance estimate of k-effective based on
! a linear combination of the collision, absorption, and tracklength
! estimates. The theory behind this can be found in M. Halperin, "Almost
! linearly-optimum combination of unbiased estimates," J. Am. Stat. Assoc., 56,
! 36-43 (1961), doi:10.1080/01621459.1961.10482088. The implementation here
@ -448,7 +449,9 @@ contains
! of keff confidence intervals in MCNP," Nucl. Technol., 111, 169-182 (1995).
!===============================================================================
subroutine calculate_combined_keff()
function openmc_get_keff(k_combined) result(err) bind(C)
real(C_DOUBLE), intent(out) :: k_combined(2)
integer(C_INT) :: err
integer :: l ! loop index
integer :: i, j, k ! indices referring to collision, absorption, or track
@ -459,9 +462,14 @@ contains
real(8) :: g ! sum of weighting factors
real(8) :: S(3) ! sums used for variance calculation
k_combined = ZERO
! Make sure we have at least four realizations. Notice that at the end,
! there is a N-3 term in a denominator.
if (n_realizations <= 3) return
if (n_realizations <= 3) then
err = -1
return
end if
! Initialize variables
n = real(n_realizations, 8)
@ -523,7 +531,6 @@ contains
! Initialize variables
g = ZERO
S = ZERO
k_combined = ZERO
do l = 1, 3
! Permutations of estimates
@ -590,8 +597,9 @@ contains
k_combined(2) = sqrt(k_combined(2))
end if
err = 0
end subroutine calculate_combined_keff
end function openmc_get_keff
!===============================================================================
! COUNT_SOURCE_FOR_UFS determines the source fraction in each UFS mesh cell and

1
src/global.F90
View file

@ -256,7 +256,6 @@ module global
real(8) :: k_col_abs = ZERO ! sum over batches of k_collision * k_absorption
real(8) :: k_col_tra = ZERO ! sum over batches of k_collision * k_tracklength
real(8) :: k_abs_tra = ZERO ! sum over batches of k_absorption * k_tracklength
real(8) :: k_combined(2) ! combined best estimate of k-effective
! Shannon entropy
logical :: entropy_on = .false.

11
src/output.F90
View file

@ -1,8 +1,10 @@
module output
use, intrinsic :: ISO_C_BINDING
use, intrinsic :: ISO_FORTRAN_ENV
use constants
use eigenvalue, only: openmc_get_keff
use endf, only: reaction_name
use error, only: fatal_error, warning
use geometry_header, only: Cell, Universe, Lattice, RectLattice, &
@ -608,6 +610,8 @@ contains
real(8) :: t_n1 ! t-value with N-1 degrees of freedom
real(8) :: t_n3 ! t-value with N-3 degrees of freedom
real(8) :: x(2) ! mean and standard deviation
real(C_DOUBLE) :: k_combined(2)
integer(C_INT) :: err
! display header block for results
call header("Results", 4)
@ -634,8 +638,11 @@ contains
write(ou,102) "k-effective (Track-length)", x(1), t_n1 * x(2)
x(:) = mean_stdev(r(:, K_ABSORPTION), n)
write(ou,102) "k-effective (Absorption)", x(1), t_n1 * x(2)
if (n > 3) write(ou,102) "Combined k-effective", k_combined(1), &
t_n3 * k_combined(2)
if (n > 3) then
err = openmc_get_keff(k_combined)
write(ou,102) "Combined k-effective", k_combined(1), &
t_n3 * k_combined(2)
end if
end if
x(:) = mean_stdev(global_tallies(:, LEAKAGE), n)
write(ou,102) "Leakage Fraction", x(1), t_n1 * x(2)

15
src/simulation.F90
View file

@ -5,9 +5,8 @@ module simulation
use cmfd_execute, only: cmfd_init_batch, execute_cmfd
use constants, only: ZERO
use eigenvalue, only: count_source_for_ufs, calculate_average_keff, &
calculate_combined_keff, calculate_generation_keff, &
shannon_entropy, synchronize_bank, keff_generation, &
k_sum
calculate_generation_keff, shannon_entropy, &
synchronize_bank, keff_generation, k_sum
#ifdef _OPENMP
use eigenvalue, only: join_bank_from_threads
#endif
@ -300,9 +299,6 @@ contains
if (run_mode == MODE_EIGENVALUE) then
! Perform CMFD calculation if on
if (cmfd_on) call execute_cmfd()
! Calculate combined estimate of k-effective
if (master) call calculate_combined_keff()
end if
! Check_triggers
@ -327,13 +323,6 @@ contains
call write_source_point()
end if
if (master .and. current_batch == n_max_batches .and. &
run_mode == MODE_EIGENVALUE) then
! Make sure combined estimate of k-effective is calculated at the last
! batch in case no state point is written
call calculate_combined_keff()
end if
end subroutine finalize_batch
!===============================================================================

8
src/state_point.F90
View file

@ -11,11 +11,12 @@ module state_point
! intervals, using the <state_point ... /> tag.
!===============================================================================
use, intrinsic :: ISO_C_BINDING, only: c_loc, c_ptr
use, intrinsic :: ISO_C_BINDING
use hdf5
use constants
use eigenvalue, only: openmc_get_keff
use endf, only: reaction_name
use error, only: fatal_error, warning
use global
@ -46,6 +47,8 @@ contains
integer(HID_T) :: cmfd_group, tallies_group, tally_group, meshes_group, &
mesh_group, filters_group, filter_group, derivs_group, &
deriv_group, runtime_group
integer(C_INT) :: err
real(C_DOUBLE) :: k_combined(2)
character(MAX_WORD_LEN), allocatable :: str_array(:)
character(MAX_FILE_LEN) :: filename
type(TallyObject), pointer :: tally
@ -117,6 +120,7 @@ contains
call write_dataset(file_id, "k_col_abs", k_col_abs)
call write_dataset(file_id, "k_col_tra", k_col_tra)
call write_dataset(file_id, "k_abs_tra", k_abs_tra)
err = openmc_get_keff(k_combined)
call write_dataset(file_id, "k_combined", k_combined)
! Write out CMFD info
@ -647,7 +651,6 @@ contains
integer(HID_T) :: cmfd_group
integer(HID_T) :: tallies_group
integer(HID_T) :: tally_group
real(8) :: real_array(3)
logical :: source_present
character(MAX_WORD_LEN) :: word
@ -727,7 +730,6 @@ contains
call read_dataset(k_col_abs, file_id, "k_col_abs")
call read_dataset(k_col_tra, file_id, "k_col_tra")
call read_dataset(k_abs_tra, file_id, "k_abs_tra")
call read_dataset(real_array(1:2), file_id, "k_combined")
! Take maximum of statepoint n_inactive and input n_inactive
n_inactive = max(n_inactive, int_array(1))

6
src/trigger.F90
View file

@ -1,10 +1,13 @@
module trigger
use, intrinsic :: ISO_C_BINDING
#ifdef MPI
use message_passing
#endif
use constants
use eigenvalue, only: openmc_get_keff
use global
use string, only: to_str
use output, only: warning, write_message
@ -101,10 +104,12 @@ contains
integer :: score_index ! scoring bin index
integer :: n_order ! loop index for moment orders
integer :: nm_order ! loop index for Ynm moment orders
integer(C_INT) :: err
real(8) :: uncertainty ! trigger uncertainty
real(8) :: std_dev = ZERO ! trigger standard deviation
real(8) :: rel_err = ZERO ! trigger relative error
real(8) :: ratio ! ratio of the uncertainty/trigger threshold
real(C_DOUBLE) :: k_combined(2)
type(TallyObject), pointer :: t ! tally pointer
type(TriggerObject), pointer :: trigger ! tally trigger
@ -119,6 +124,7 @@ contains
! Check eigenvalue trigger
if (run_mode == MODE_EIGENVALUE) then
if (keff_trigger % trigger_type /= 0) then
err = openmc_get_keff(k_combined)
select case (keff_trigger % trigger_type)
case(VARIANCE)
uncertainty = k_combined(2) ** 2