crp/OpenMC
1
0
Fork 0
mirror of https://github.com/openmc-dev/openmc.git synced 2026-07-27 05:35:49 -04:00
Code Issues Projects Releases Packages Wiki Activity Actions

Add/use Fortran bindings to ThermalScattering C++ class

Browse source
This commit is contained in:
Paul Romano 2018-08-15 12:59:20 -05:00
parent c4680b43b7
commit b7f1f7cf4b
6 changed files with 188 additions and 674 deletions
Download patch file Download diff file Expand all files Collapse all files

4
src/material_header.F90
View file

@ -172,7 +172,7 @@ contains
found = .false.
associate (sab => sab_tables(this % i_sab_tables(k)))
FIND_NUCLIDE: do j = 1, size(this % nuclide)
if (any(sab % nuclides == nuclides(this % nuclide(j)) % name)) then
if (sab % has_nuclide(nuclides(this % nuclide(j)) % name)) then
call i_sab_tables % push_back(this % i_sab_tables(k))
call i_sab_nuclides % push_back(j)
call sab_fracs % push_back(this % sab_fracs(k))
@ -326,7 +326,7 @@ contains
! If particle energy is greater than the highest energy for the
! S(a,b) table, then don't use the S(a,b) table
if (p % E > sab_tables(i_sab) % data(1) % threshold_inelastic) then
if (p % E > sab_tables(i_sab) % threshold()) then
i_sab = 0
end if

8
src/nuclide_header.F90
View file

@ -198,7 +198,7 @@ module nuclide_header
type(DictCharInt) :: nuclide_dict
! Cross section caches
type(NuclideMicroXS), allocatable :: micro_xs(:) ! Cache for each nuclide
type(NuclideMicroXS), allocatable, target :: micro_xs(:) ! Cache for each nuclide
type(MaterialMacroXS) :: material_xs ! Cache for current material
!$omp threadprivate(micro_xs, material_xs)
@ -1096,9 +1096,9 @@ contains
real(8), intent(in) :: sab_frac ! fraction of atoms affected by S(a,b)
type(NuclideMicroXS), intent(inout) :: micro_xs ! Cross section cache
integer :: i_temp ! temperature index
real(8) :: inelastic ! S(a,b) inelastic cross section
real(8) :: elastic ! S(a,b) elastic cross section
integer(C_INT) :: i_temp ! temperature index
real(C_DOUBLE) :: inelastic ! S(a,b) inelastic cross section
real(C_DOUBLE) :: elastic ! S(a,b) elastic cross section
! Set flag that S(a,b) treatment should be used for scattering
micro_xs % index_sab = i_sab

261
src/physics.F90
View file

@ -884,263 +884,16 @@ contains
real(8), intent(inout) :: uvw(3) ! directional cosines
real(8), intent(out) :: mu ! scattering cosine
integer :: i ! incoming energy bin
integer :: j ! outgoing energy bin
integer :: k ! outgoing cosine bin
integer :: i_temp ! temperature index
integer :: n_energy_out ! number of outgoing energy bins
real(8) :: f ! interpolation factor
real(8) :: r ! used for skewed sampling & continuous
real(8) :: E_ij ! outgoing energy j for E_in(i)
real(8) :: E_i1j ! outgoing energy j for E_in(i+1)
real(8) :: mu_ijk ! outgoing cosine k for E_in(i) and E_out(j)
real(8) :: mu_i1jk ! outgoing cosine k for E_in(i+1) and E_out(j)
real(8) :: prob ! probability for sampling Bragg edge
! Following are needed only for SAB_SECONDARY_CONT scattering
integer :: l ! sampled incoming E bin (is i or i + 1)
real(8) :: E_i_1, E_i_J ! endpoints on outgoing grid i
real(8) :: E_i1_1, E_i1_J ! endpoints on outgoing grid i+1
real(8) :: E_1, E_J ! endpoints interpolated between i and i+1
real(8) :: E_l_j, E_l_j1 ! adjacent E on outgoing grid l
real(8) :: p_l_j, p_l_j1 ! adjacent p on outgoing grid l
real(8) :: c_j, c_j1 ! cumulative probability
real(8) :: frac ! interpolation factor on outgoing energy
real(8) :: r1 ! RNG for outgoing energy
real(8) :: mu_left, mu_right ! adjacent mu values
real(C_DOUBLE) :: E_out
type(C_PTR) :: ptr
i_temp = micro_xs(i_nuclide) % index_temp_sab
! Sample from C++ side
ptr = C_LOC(micro_xs(i_nuclide))
call sab_tables(i_sab) % sample(ptr, E, E_out, mu)
! Get pointer to S(a,b) table
associate (sab => sab_tables(i_sab) % data(i_temp))
! Determine whether inelastic or elastic scattering will occur
if (prn() < micro_xs(i_nuclide) % thermal_elastic / &
micro_xs(i_nuclide) % thermal) then
! elastic scattering
! Get index and interpolation factor for elastic grid
if (E < sab % elastic_e_in(1)) then
i = 1
f = ZERO
else
i = binary_search(sab % elastic_e_in, sab % n_elastic_e_in, E)
f = (E - sab%elastic_e_in(i)) / &
(sab%elastic_e_in(i+1) - sab%elastic_e_in(i))
end if
! Select treatment based on elastic mode
if (sab % elastic_mode == SAB_ELASTIC_DISCRETE) then
! With this treatment, we interpolate between two discrete cosines
! corresponding to neighboring incoming energies. This is used for
! data derived in the incoherent approximation
! Sample outgoing cosine bin
k = 1 + int(prn() * sab % n_elastic_mu)
! Determine outgoing cosine corresponding to E_in(i) and E_in(i+1)
mu_ijk = sab % elastic_mu(k,i)
mu_i1jk = sab % elastic_mu(k,i+1)
! Cosine of angle between incoming and outgoing neutron
mu = (1 - f)*mu_ijk + f*mu_i1jk
elseif (sab % elastic_mode == SAB_ELASTIC_EXACT) then
! This treatment is used for data derived in the coherent
! approximation, i.e. for crystalline structures that have Bragg
! edges.
! Sample a Bragg edge between 1 and i
prob = prn() * sab % elastic_P(i+1)
if (prob < sab % elastic_P(1)) then
k = 1
else
k = binary_search(sab % elastic_P(1:i+1), i+1, prob)
end if
! Characteristic scattering cosine for this Bragg edge
mu = ONE - TWO*sab % elastic_e_in(k) / E
end if
! Outgoing energy is same as incoming energy -- no need to do anything
else
! Perform inelastic calculations
! Get index and interpolation factor for inelastic grid
if (E < sab % inelastic_e_in(1)) then
i = 1
f = ZERO
else
i = binary_search(sab % inelastic_e_in, sab % n_inelastic_e_in, E)
f = (E - sab%inelastic_e_in(i)) / &
(sab%inelastic_e_in(i+1) - sab%inelastic_e_in(i))
end if
! Now that we have an incoming energy bin, we need to determine the
! outgoing energy bin. This will depend on the "secondary energy
! mode". If the mode is 0, then the outgoing energy bin is chosen from a
! set of equally-likely bins. If the mode is 1, then the first
! two and last two bins are skewed to have lower probabilities than the
! other bins (0.1 for the first and last bins and 0.4 for the second and
! second to last bins, relative to a normal bin probability of 1).
! Finally, if the mode is 2, then a continuous distribution (with
! accompanying PDF and CDF is utilized)
if ((sab_tables(i_sab) % secondary_mode == SAB_SECONDARY_EQUAL) .or. &
(sab_tables(i_sab) % secondary_mode == SAB_SECONDARY_SKEWED)) then
if (sab_tables(i_sab) % secondary_mode == SAB_SECONDARY_EQUAL) then
! All bins equally likely
j = 1 + int(prn() * sab % n_inelastic_e_out)
elseif (sab_tables(i_sab) % secondary_mode == SAB_SECONDARY_SKEWED) then
! Distribution skewed away from edge points
! Determine number of outgoing energy and angle bins
n_energy_out = sab % n_inelastic_e_out
r = prn() * (n_energy_out - 3)
if (r > ONE) then
! equally likely N-4 middle bins
j = int(r) + 2
elseif (r > 0.6_8) then
! second to last bin has relative probability of 0.4
j = n_energy_out - 1
elseif (r > HALF) then
! last bin has relative probability of 0.1
j = n_energy_out
elseif (r > 0.1_8) then
! second bin has relative probability of 0.4
j = 2
else
! first bin has relative probability of 0.1
j = 1
end if
end if
! Determine outgoing energy corresponding to E_in(i) and E_in(i+1)
E_ij = sab % inelastic_e_out(j,i)
E_i1j = sab % inelastic_e_out(j,i+1)
! Outgoing energy
E = (1 - f)*E_ij + f*E_i1j
! Sample outgoing cosine bin
k = 1 + int(prn() * sab % n_inelastic_mu)
! Determine outgoing cosine corresponding to E_in(i) and E_in(i+1)
mu_ijk = sab % inelastic_mu(k,j,i)
mu_i1jk = sab % inelastic_mu(k,j,i+1)
! Cosine of angle between incoming and outgoing neutron
mu = (1 - f)*mu_ijk + f*mu_i1jk
else if (sab_tables(i_sab) % secondary_mode == SAB_SECONDARY_CONT) then
! Continuous secondary energy - this is to be similar to
! Law 61 interpolation on outgoing energy
! Sample between ith and (i+1)th bin
r = prn()
if (f > r) then
l = i + 1
else
l = i
end if
! Determine endpoints on grid i
n_energy_out = sab % inelastic_data(i) % n_e_out
E_i_1 = sab % inelastic_data(i) % e_out(1)
E_i_J = sab % inelastic_data(i) % e_out(n_energy_out)
! Determine endpoints on grid i + 1
n_energy_out = sab % inelastic_data(i + 1) % n_e_out
E_i1_1 = sab % inelastic_data(i + 1) % e_out(1)
E_i1_J = sab % inelastic_data(i + 1) % e_out(n_energy_out)
E_1 = E_i_1 + f * (E_i1_1 - E_i_1)
E_J = E_i_J + f * (E_i1_J - E_i_J)
! Determine outgoing energy bin
! (First reset n_energy_out to the right value)
n_energy_out = sab % inelastic_data(l) % n_e_out
r1 = prn()
c_j = sab % inelastic_data(l) % e_out_cdf(1)
do j = 1, n_energy_out - 1
c_j1 = sab % inelastic_data(l) % e_out_cdf(j + 1)
if (r1 < c_j1) exit
c_j = c_j1
end do
! check to make sure k is <= n_energy_out - 1
j = min(j, n_energy_out - 1)
! Get the data to interpolate between
E_l_j = sab % inelastic_data(l) % e_out(j)
p_l_j = sab % inelastic_data(l) % e_out_pdf(j)
! Next part assumes linear-linear interpolation in standard
E_l_j1 = sab % inelastic_data(l) % e_out(j + 1)
p_l_j1 = sab % inelastic_data(l) % e_out_pdf(j + 1)
! Find secondary energy (variable E)
frac = (p_l_j1 - p_l_j) / (E_l_j1 - E_l_j)
if (frac == ZERO) then
E = E_l_j + (r1 - c_j) / p_l_j
else
E = E_l_j + (sqrt(max(ZERO, p_l_j * p_l_j + &
TWO * frac * (r1 - c_j))) - p_l_j) / frac
end if
! Now interpolate between incident energy bins i and i + 1
if (l == i) then
E = E_1 + (E - E_i_1) * (E_J - E_1) / (E_i_J - E_i_1)
else
E = E_1 + (E - E_i1_1) * (E_J - E_1) / (E_i1_J - E_i1_1)
end if
! Sample outgoing cosine bin
k = 1 + int(prn() * sab % n_inelastic_mu)
! Rather than use the sampled discrete mu directly, it is smeared over
! a bin of width min(mu[k] - mu[k-1], mu[k+1] - mu[k]) centered on the
! discrete mu value itself.
associate (mu_l => sab % inelastic_data(l) % mu)
f = (r1 - c_j)/(c_j1 - c_j)
! Determine (k-1)th mu value
if (k == 1) then
mu_left = -ONE
else
mu_left = mu_l(k-1, j) + f*(mu_l(k-1, j+1) - mu_l(k-1,j))
end if
! Determine kth mu value
mu = mu_l(k, j) + f*(mu_l(k, j+1) - mu_l(k, j))
! Determine (k+1)th mu value
if (k == sab % n_inelastic_mu) then
mu_right = ONE - mu
else
mu_right = mu_l(k+1, j) + f*(mu_l(k+1, j+1) - mu_l(k+1,j)) - mu
end if
end associate
! Smear angle
mu = mu + min(mu - mu_left, mu_right - mu)*(prn() - HALF)
end if ! (inelastic secondary energy treatment)
end if ! (elastic or inelastic)
end associate
! Because of floating-point roundoff, it may be possible for mu to be
! outside of the range [-1,1). In these cases, we just set mu to exactly
! -1 or 1
if (abs(mu) > ONE) mu = sign(ONE,mu)
! change direction of particle
! Set energy to outgoing, change direction of particle
E = E_out
uvw = rotate_angle(uvw, mu)
end subroutine sab_scatter
!===============================================================================

530
src/sab_header.F90
View file

@ -1,467 +1,169 @@
module sab_header
use, intrinsic :: ISO_C_BINDING
use, intrinsic :: ISO_FORTRAN_ENV
use algorithm, only: find, sort, binary_search
use constants
use dict_header, only: DictIntInt, DictCharInt
use distribution_univariate, only: Tabular
use error, only: warning, fatal_error
use dict_header, only: DictCharInt
use hdf5_interface
use random_lcg, only: prn
use secondary_correlated, only: CorrelatedAngleEnergy
use settings
use stl_vector, only: VectorInt, VectorReal
use string, only: to_str, str_to_int
use stl_vector, only: VectorReal
use string, only: to_c_string
implicit none
private
!===============================================================================
! DISTENERGYSAB contains the secondary energy/angle distributions for inelastic
! thermal scattering collisions which utilize a continuous secondary energy
! representation.
!===============================================================================
type DistEnergySab
integer :: n_e_out
real(8), allocatable :: e_out(:)
real(8), allocatable :: e_out_pdf(:)
real(8), allocatable :: e_out_cdf(:)
real(8), allocatable :: mu(:,:)
end type DistEnergySab
public :: free_memory_sab
!===============================================================================
! SALPHABETA contains S(a,b) data for thermal neutron scattering, typically off
! of light isotopes such as water, graphite, Be, etc
!===============================================================================
type SabData
! threshold for S(a,b) treatment (usually ~4 eV)
real(8) :: threshold_inelastic
real(8) :: threshold_elastic = ZERO
! Inelastic scattering data
integer :: n_inelastic_e_in ! # of incoming E for inelastic
integer :: n_inelastic_e_out ! # of outgoing E for inelastic
integer :: n_inelastic_mu ! # of outgoing angles for inelastic
real(8), allocatable :: inelastic_e_in(:)
real(8), allocatable :: inelastic_sigma(:)
! The following are used only if secondary_mode is 0 or 1
real(8), allocatable :: inelastic_e_out(:,:)
real(8), allocatable :: inelastic_mu(:,:,:)
! The following is used only if secondary_mode is 3
! The different implementation is necessary because the continuous
! representation has a variable number of outgoing energy points for each
! incoming energy
type(DistEnergySab), allocatable :: inelastic_data(:) ! One for each Ein
! Elastic scattering data
integer :: elastic_mode ! elastic mode (discrete/exact)
integer :: n_elastic_e_in ! # of incoming E for elastic
integer :: n_elastic_mu ! # of outgoing angles for elastic
real(8), allocatable :: elastic_e_in(:)
real(8), allocatable :: elastic_P(:)
real(8), allocatable :: elastic_mu(:,:)
end type SabData
type SAlphaBeta
character(150) :: name ! name of table, e.g. lwtr.10t
real(8) :: awr ! weight of nucleus in neutron masses
real(8), allocatable :: kTs(:) ! temperatures in eV (k*T)
character(10), allocatable :: nuclides(:) ! List of valid nuclides
integer :: secondary_mode ! secondary mode (equal/skewed/continuous)
! cross sections and distributions at each temperature
type(SabData), allocatable :: data(:)
type, public :: SAlphaBeta
type(C_PTR) :: ptr
contains
procedure :: from_hdf5 => salphabeta_from_hdf5
procedure :: calculate_xs => sab_calculate_xs
procedure :: from_hdf5
procedure :: calculate_xs
procedure :: free
procedure :: has_nuclide
procedure :: sample
procedure :: threshold
end type SAlphaBeta
! S(a,b) tables
type(SAlphaBeta), allocatable, target :: sab_tables(:)
integer(C_INT), bind(C) :: n_sab_tables
type(DictCharInt) :: sab_dict
type(SAlphaBeta), public, allocatable, target :: sab_tables(:)
integer(C_INT), public, bind(C) :: n_sab_tables
type(DictCharInt), public :: sab_dict
interface
function sab_from_hdf5(group_id, temperature, n, method, &
tolerance, minmax) result(ptr) bind(C)
import HID_T, C_DOUBLE, C_INT, C_PTR
integer(HID_T), value :: group_id
real(C_DOUBLE), intent(in) :: temperature
integer(C_INT), value :: n
integer(C_INT), value :: method
real(C_DOUBLE), value :: tolerance
real(C_DOUBLE), intent(in) :: minmax(2)
type(C_PTR) :: ptr
end function
subroutine sab_calculate_xs(ptr, E, sqrtkT, i_temp, elastic, &
inelastic) bind(C)
import C_PTR, C_DOUBLE, C_INT
type(C_PTR), value :: ptr
real(C_DOUBLE), value :: E
real(C_DOUBLE), value :: sqrtkT
integer(C_INT), intent(out) :: i_temp
real(C_DOUBLE), intent(out) :: elastic
real(C_DOUBLE), intent(out) :: inelastic
end subroutine
subroutine sab_free(ptr) bind(C)
import C_PTR
type(C_PTR), value :: ptr
end subroutine
function sab_has_nuclide(ptr, name) result(val) bind(C)
import C_PTR, C_CHAR, C_BOOL
type(C_PTR), value :: ptr
character(kind=C_CHAR), intent(in) :: name(*)
logical(C_BOOL) :: val
end function
subroutine sab_sample(ptr, micro_xs, E_in, E_out, mu) bind(C)
import C_PTR, C_INT, C_DOUBLE
type(C_PTR), value :: ptr
type(C_PTR), value :: micro_xs
real(C_DOUBLE), value :: E_in
real(C_DOUBLE), intent(out) :: E_out
real(C_DOUBLE), intent(out) :: mu
end subroutine
function sab_threshold(ptr) result(threshold) bind(C)
import C_PTR, C_double
type(C_PTR), value :: ptr
real(C_DOUBLE) :: threshold
end function
end interface
contains
subroutine salphabeta_from_hdf5(this, group_id, temperature, method, &
subroutine from_hdf5(this, group_id, temperature, method, &
tolerance, minmax)
class(SAlphaBeta), intent(inout) :: this
integer(HID_T), intent(in) :: group_id
type(VectorReal), intent(in) :: temperature ! list of temperatures
integer, intent(in) :: method
real(8), intent(in) :: tolerance
real(8), intent(in) :: minmax(2)
integer(C_INT), intent(in) :: method
real(C_DOUBLE), intent(in) :: tolerance
real(C_DOUBLE), intent(in) :: minmax(2)
integer :: i, j
integer :: t
integer :: n_energy, n_energy_out, n_mu
integer :: i_closest
integer :: n_temperature
integer(HID_T) :: T_group
integer(HID_T) :: elastic_group
integer(HID_T) :: inelastic_group
integer(HID_T) :: dset_id
integer(HID_T) :: kT_group
integer(HSIZE_T) :: dims2(2)
integer(HSIZE_T) :: dims3(3)
real(8), allocatable :: temp(:,:)
character(20) :: type
type(CorrelatedAngleEnergy) :: correlated_dist
real(C_DOUBLE) :: dummy
integer(C_INT) :: n
character(MAX_WORD_LEN) :: temp_str
character(MAX_WORD_LEN), allocatable :: dset_names(:)
real(8), allocatable :: temps_available(:) ! temperatures available
real(8) :: temp_desired
real(8) :: temp_actual
type(VectorInt) :: temps_to_read
! Get name of table from group
this % name = get_name(group_id)
! Get rid of leading '/'
this % name = trim(this % name(2:))
call read_attribute(this % awr, group_id, 'atomic_weight_ratio')
call read_attribute(this % nuclides, group_id, 'nuclides')
call read_attribute(type, group_id, 'secondary_mode')
select case (type)
case ('equal')
this % secondary_mode = SAB_SECONDARY_EQUAL
case ('skewed')
this % secondary_mode = SAB_SECONDARY_SKEWED
case ('continuous')
this % secondary_mode = SAB_SECONDARY_CONT
end select
! Read temperatures
kT_group = open_group(group_id, 'kTs')
! Determine temperatures available
call get_datasets(kT_group, dset_names)
allocate(temps_available(size(dset_names)))
do i = 1, size(dset_names)
! Read temperature value
call read_dataset(temps_available(i), kT_group, trim(dset_names(i)))
temps_available(i) = temps_available(i) / K_BOLTZMANN
end do
call sort(temps_available)
! Determine actual temperatures to read -- start by checking whether a
! temperature range was given, in which case all temperatures in the range
! are loaded irrespective of what temperatures actually appear in the model
if (minmax(2) > ZERO) then
do i = 1, size(temps_available)
temp_actual = temps_available(i)
if (minmax(1) <= temp_actual .and. temp_actual <= minmax(2)) then
call temps_to_read % push_back(nint(temp_actual))
end if
end do
n = temperature % size()
if (n > 0) then
this % ptr = sab_from_hdf5(group_id, temperature % data(1), n, method, tolerance, minmax)
else
! In this case, temperatures % data(1) doesn't exist, so we just pass a
! dummy value
this % ptr = sab_from_hdf5(group_id, dummy, n, method, tolerance, minmax)
end if
select case (method)
case (TEMPERATURE_NEAREST)
! Determine actual temperatures to read
do i = 1, temperature % size()
temp_desired = temperature % data(i)
i_closest = minloc(abs(temps_available - temp_desired), dim=1)
temp_actual = temps_available(i_closest)
if (abs(temp_actual - temp_desired) < tolerance) then
if (find(temps_to_read, nint(temp_actual)) == -1) then
call temps_to_read % push_back(nint(temp_actual))
end if
else
call fatal_error("Nuclear data library does not contain cross sections &
&for " // trim(this % name) // " at or near " // &
trim(to_str(nint(temp_desired))) // " K.")
end if
end do
case (TEMPERATURE_INTERPOLATION)
! If temperature interpolation or multipole is selected, get a list of
! bounding temperatures for each actual temperature present in the model
TEMP_LOOP: do i = 1, temperature % size()
temp_desired = temperature % data(i)
do j = 1, size(temps_available) - 1
if (temps_available(j) <= temp_desired .and. &
temp_desired < temps_available(j + 1)) then
if (find(temps_to_read, nint(temps_available(j))) == -1) then
call temps_to_read % push_back(nint(temps_available(j)))
end if
if (find(temps_to_read, nint(temps_available(j + 1))) == -1) then
call temps_to_read % push_back(nint(temps_available(j + 1)))
end if
cycle TEMP_LOOP
end if
end do
call fatal_error("Nuclear data library does not contain cross sections &
&for " // trim(this % name) // " at temperatures that bound " // &
trim(to_str(nint(temp_desired))) // " K.")
end do TEMP_LOOP
end select
! Sort temperatures to read
call sort(temps_to_read)
n_temperature = temps_to_read % size()
allocate(this % kTs(n_temperature))
allocate(this % data(n_temperature))
do t = 1, n_temperature
! Get temperature as a string
temp_str = trim(to_str(temps_to_read % data(t))) // "K"
! Read exact temperature value
call read_dataset(this % kTs(t), kT_group, temp_str)
! Open group for temperature i
T_group = open_group(group_id, temp_str)
! Coherent elastic data
if (object_exists(T_group, 'elastic')) then
! Read cross section data
elastic_group = open_group(T_group, 'elastic')
dset_id = open_dataset(elastic_group, 'xs')
call read_attribute(type, dset_id, 'type')
call get_shape(dset_id, dims2)
allocate(temp(dims2(1), dims2(2)))
call read_dataset(temp, dset_id)
call close_dataset(dset_id)
! Set cross section data and type
this % data(t) % n_elastic_e_in = int(dims2(1), 4)
allocate(this % data(t) % elastic_e_in(this % data(t) % n_elastic_e_in))
allocate(this % data(t) % elastic_P(this % data(t) % n_elastic_e_in))
this % data(t) % elastic_e_in(:) = temp(:, 1)
this % data(t) % elastic_P(:) = temp(:, 2)
select case (type)
case ('tab1')
this % data(t) % elastic_mode = SAB_ELASTIC_DISCRETE
case ('bragg')
this % data(t) % elastic_mode = SAB_ELASTIC_EXACT
end select
deallocate(temp)
! Set elastic threshold
this % data(t) % threshold_elastic = this % data(t) % elastic_e_in(&
this % data(t) % n_elastic_e_in)
! Read angle distribution
if (this % data(t) % elastic_mode /= SAB_ELASTIC_EXACT) then
dset_id = open_dataset(elastic_group, 'mu_out')
call get_shape(dset_id, dims2)
this % data(t) % n_elastic_mu = int(dims2(1), 4)
allocate(this % data(t) % elastic_mu(dims2(1), dims2(2)))
call read_dataset(this % data(t) % elastic_mu, dset_id)
call close_dataset(dset_id)
end if
call close_group(elastic_group)
end if
! Inelastic data
if (object_exists(T_group, 'inelastic')) then
! Read type of inelastic data
inelastic_group = open_group(T_group, 'inelastic')
! Read cross section data
dset_id = open_dataset(inelastic_group, 'xs')
call get_shape(dset_id, dims2)
allocate(temp(dims2(1), dims2(2)))
call read_dataset(temp, dset_id)
call close_dataset(dset_id)
! Set cross section data
this % data(t) % n_inelastic_e_in = int(dims2(1), 4)
allocate(this % data(t) % inelastic_e_in(this % data(t) % n_inelastic_e_in))
allocate(this % data(t) % inelastic_sigma(this % data(t) % n_inelastic_e_in))
this % data(t) % inelastic_e_in(:) = temp(:, 1)
this % data(t) % inelastic_sigma(:) = temp(:, 2)
deallocate(temp)
! Set inelastic threshold
this % data(t) % threshold_inelastic = this % data(t) % inelastic_e_in(&
this % data(t) % n_inelastic_e_in)
if (this % secondary_mode /= SAB_SECONDARY_CONT) then
! Read energy distribution
dset_id = open_dataset(inelastic_group, 'energy_out')
call get_shape(dset_id, dims2)
this % data(t) % n_inelastic_e_out = int(dims2(1), 4)
allocate(this % data(t) % inelastic_e_out(dims2(1), dims2(2)))
call read_dataset(this % data(t) % inelastic_e_out, dset_id)
call close_dataset(dset_id)
! Read angle distribution
dset_id = open_dataset(inelastic_group, 'mu_out')
call get_shape(dset_id, dims3)
this % data(t) % n_inelastic_mu = int(dims3(1), 4)
allocate(this % data(t) % inelastic_mu(dims3(1), dims3(2), dims3(3)))
call read_dataset(this % data(t) % inelastic_mu, dset_id)
call close_dataset(dset_id)
else
! Read correlated angle-energy distribution
call correlated_dist % from_hdf5(inelastic_group)
! Convert to S(a,b) native format
n_energy = size(correlated_dist % energy)
allocate(this % data(t) % inelastic_data(n_energy))
do i = 1, n_energy
associate (edist => correlated_dist % distribution(i))
! Get number of outgoing energies for incoming energy i
n_energy_out = size(edist % e_out)
this % data(t) % inelastic_data(i) % n_e_out = n_energy_out
allocate(this % data(t) % inelastic_data(i) % e_out(n_energy_out))
allocate(this % data(t) % inelastic_data(i) % e_out_pdf(n_energy_out))
allocate(this % data(t) % inelastic_data(i) % e_out_cdf(n_energy_out))
! Copy outgoing energy distribution
this % data(t) % inelastic_data(i) % e_out(:) = edist % e_out
this % data(t) % inelastic_data(i) % e_out_pdf(:) = edist % p
this % data(t) % inelastic_data(i) % e_out_cdf(:) = edist % c
do j = 1, n_energy_out
select type (adist => edist % angle(j) % obj)
type is (Tabular)
! On first pass, allocate space for angles
if (j == 1) then
n_mu = size(adist % x)
this % data(t) % n_inelastic_mu = n_mu
allocate(this % data(t) % inelastic_data(i) % mu(&
n_mu, n_energy_out))
end if
! Copy outgoing angles
this % data(t) % inelastic_data(i) % mu(:, j) = adist % x
end select
end do
end associate
end do
! Clear data on correlated angle-energy object
deallocate(correlated_dist % breakpoints)
deallocate(correlated_dist % interpolation)
deallocate(correlated_dist % energy)
deallocate(correlated_dist % distribution)
end if
call close_group(inelastic_group)
end if
call close_group(T_group)
end do
call close_group(kT_group)
end subroutine salphabeta_from_hdf5
end subroutine from_hdf5
!===============================================================================
! SAB_CALCULATE_XS determines the elastic and inelastic scattering
! cross-sections in the thermal energy range.
!===============================================================================
subroutine sab_calculate_xs(this, E, sqrtkT, i_temp, elastic, inelastic)
subroutine calculate_xs(this, E, sqrtkT, i_temp, elastic, inelastic)
class(SAlphaBeta), intent(in) :: this ! S(a,b) object
real(8), intent(in) :: E ! energy
real(8), intent(in) :: sqrtkT ! temperature
real(C_DOUBLE), intent(in) :: E ! energy
real(C_DOUBLE), intent(in) :: sqrtkT ! temperature
integer, intent(out) :: i_temp ! index in the S(a,b)'s temperature
real(8), intent(out) :: elastic ! thermal elastic cross section
real(8), intent(out) :: inelastic ! thermal inelastic cross section
real(C_DOUBLE), intent(out) :: elastic ! thermal elastic cross section
real(C_DOUBLE), intent(out) :: inelastic ! thermal inelastic cross section
integer :: i_grid ! index on S(a,b) energy grid
real(8) :: f ! interp factor on S(a,b) energy grid
real(8) :: kT
call sab_calculate_xs(this % ptr, E, sqrtkT, i_temp, elastic, inelastic)
end subroutine
! Determine temperature for S(a,b) table
kT = sqrtkT**2
if (temperature_method == TEMPERATURE_NEAREST) then
! If using nearest temperature, do linear search on temperature
do i_temp = 1, size(this % kTs)
if (abs(this % kTs(i_temp) - kT) < &
K_BOLTZMANN*temperature_tolerance) exit
end do
else
! Find temperatures that bound the actual temperature
do i_temp = 1, size(this % kTs) - 1
if (this % kTs(i_temp) <= kT .and. &
kT < this % kTs(i_temp + 1)) exit
end do
subroutine free(this)
class(SAlphaBeta), intent(inout) :: this
call sab_free(this % ptr)
end subroutine
! Randomly sample between temperature i and i+1
f = (kT - this % kTs(i_temp)) / &
(this % kTs(i_temp + 1) - this % kTs(i_temp))
if (f > prn()) i_temp = i_temp + 1
end if
function has_nuclide(this, name) result(val)
class(SAlphaBeta), intent(in) :: this
character(len=*), intent(in) :: name
logical(C_BOOL) :: val
val = sab_has_nuclide(this % ptr, to_c_string(name))
end function
! Get pointer to S(a,b) table
associate (sab => this % data(i_temp))
subroutine sample(this, micro_xs, E_in, E_out, mu)
class(SAlphaBeta), intent(in) :: this
type(C_PTR), value :: micro_xs
real(C_DOUBLE), value :: E_in
real(C_DOUBLE), intent(out) :: E_out
real(C_DOUBLE), intent(out) :: mu
! Get index and interpolation factor for inelastic grid
if (E < sab % inelastic_e_in(1)) then
i_grid = 1
f = ZERO
else
i_grid = binary_search(sab % inelastic_e_in, sab % n_inelastic_e_in, E)
f = (E - sab%inelastic_e_in(i_grid)) / &
(sab%inelastic_e_in(i_grid+1) - sab%inelastic_e_in(i_grid))
end if
! Calculate S(a,b) inelastic scattering cross section
inelastic = (ONE - f) * sab % inelastic_sigma(i_grid) + &
f * sab % inelastic_sigma(i_grid + 1)
! Check for elastic data
if (E < sab % threshold_elastic) then
! Determine whether elastic scattering is given in the coherent or
! incoherent approximation. For coherent, the cross section is
! represented as P/E whereas for incoherent, it is simply P
if (sab % elastic_mode == SAB_ELASTIC_EXACT) then
if (E < sab % elastic_e_in(1)) then
! If energy is below that of the lowest Bragg peak, the elastic
! cross section will be zero
elastic = ZERO
else
i_grid = binary_search(sab % elastic_e_in, &
sab % n_elastic_e_in, E)
elastic = sab % elastic_P(i_grid) / E
end if
else
! Determine index on elastic energy grid
if (E < sab % elastic_e_in(1)) then
i_grid = 1
else
i_grid = binary_search(sab % elastic_e_in, &
sab % n_elastic_e_in, E)
end if
! Get interpolation factor for elastic grid
f = (E - sab%elastic_e_in(i_grid))/(sab%elastic_e_in(i_grid+1) - &
sab%elastic_e_in(i_grid))
! Calculate S(a,b) elastic scattering cross section
elastic = (ONE - f) * sab % elastic_P(i_grid) + &
f * sab % elastic_P(i_grid + 1)
end if
else
! No elastic data
elastic = ZERO
end if
end associate
end subroutine sab_calculate_xs
call sab_sample(this % ptr, micro_xs, E_in, E_out, mu)
end subroutine
function threshold(this)
class(SAlphaBeta), intent(in) :: this
real(C_DOUBLE) :: threshold
threshold = sab_threshold(this % ptr)
end function
!===============================================================================
! FREE_MEMORY_SAB deallocates global arrays defined in this module
!===============================================================================
subroutine free_memory_sab()
integer :: i
n_sab_tables = 0
do i = 1, size(sab_tables)
call sab_tables(i) % free()
end do
if (allocated(sab_tables)) deallocate(sab_tables)
call sab_dict % clear()
end subroutine free_memory_sab

37
src/thermal.cpp
View file

@ -577,4 +577,41 @@ ThermalData::sample(const NuclideMicroXS* micro_xs, double E,
}
//==============================================================================
// Fortran compatibility functions
//==============================================================================
ThermalScattering*
sab_from_hdf5(hid_t group, const double* temperature, int n,
int method, double tolerance, const double* minmax)
{
// Convert temperatures to a vector
std::vector<double> T {temperature, temperature + n};
// Create new object and return it
return new ThermalScattering{group, T, method, tolerance, minmax};
}
void sab_calculate_xs(ThermalScattering* data, double E, double sqrtkT,
int* i_temp, double* elastic, double* inelastic)
{
data->calculate_xs(E, sqrtkT, i_temp, elastic, inelastic);
}
void sab_free(ThermalScattering* data) { delete data; }
bool sab_has_nuclide(ThermalScattering* data, const char* name)
{
return data->has_nuclide(name);
}
void sab_sample(ThermalScattering* data, const NuclideMicroXS* micro_xs,
double E_in, double* E_out, double* mu)
{
int i_temp = micro_xs->index_temp_sab;
data->data_[i_temp - 1].sample(micro_xs, E_in, E_out, mu);
}
double sab_threshold(ThermalScattering* data) { return data->threshold(); }
} // namespace openmc

22
src/thermal.h
View file

@ -111,6 +111,12 @@ public:
//! \return Whether table applies to the nuclide
bool has_nuclide(const char* name) const;
// Sample an outgoing energy and angle
void sample(const NuclideMicroXS* micro_xs, double E_in,
double* E_out, double* mu);
double threshold() const { return data_[0].threshold_inelastic_; }
std::string name_; //!< name of table, e.g. "c_H_in_H2O"
double awr_; //!< weight of nucleus in neutron masses
std::vector<double> kTs_; //!< temperatures in [eV] (k*T)
@ -120,6 +126,22 @@ public:
std::vector<ThermalData> data_;
};
//==============================================================================
// Fortran compatibility functions
//==============================================================================
extern "C" {
ThermalScattering* sab_from_hdf5(hid_t group, const double* temperature,
int n, int method, double tolerance, const double* minmax);
void sab_calculate_xs(ThermalScattering* data, double E, double sqrtkT,
int* i_temp, double* elastic, double* inelastic);
void sab_free(ThermalScattering* data);
bool sab_has_nuclide(ThermalScattering* data, const char* name);
void sab_sample(ThermalScattering* data, const NuclideMicroXS* micro_xs,
double E_in, double* E_out, double* mu);
double sab_threshold(ThermalScattering* data);
}
} // namespace openmc
#endif // OPENMC_THERMAL_SCATTERING_H