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

Use Reaction type hierarchy on C++ side

Browse source
This commit is contained in:
Paul Romano 2018-07-11 13:21:05 -05:00
parent d4775552e0
commit 68c1358565
10 changed files with 502 additions and 192 deletions
Download patch file Download diff file Expand all files Collapse all files

5
src/endf.cpp
View file

@ -28,6 +28,11 @@ Interpolation int2interp(int i)
}
}
bool is_fission(int mt)
{
return mt == 18 || mt == 19 || mt == 20 || mt == 21 || mt == 38;
}
//==============================================================================
// Polynomial implementation
//==============================================================================

1
src/endf.h
View file

@ -9,6 +9,7 @@
namespace openmc {
Interpolation int2interp(int i);
bool is_fission(int MT);
class Function1D {
public:

146
src/nuclide_header.F90
View file

@ -642,29 +642,33 @@ contains
if (rx % MT >= N_2N0 .and. rx % MT <= N_2NC .and. find(MTs, N_2N) /= -1) cycle
do t = 1, n_temperature
j = rx % xs(t) % threshold
n = size(rx % xs(t) % value)
j = rx % xs_threshold(t)
n = rx % xs_size(t)
! Add contribution to total cross section
this % xs(t) % value(XS_TOTAL,j:j+n-1) = this % xs(t) % &
value(XS_TOTAL,j:j+n-1) + rx % xs(t) % value
do k = j, j + n - 1
this % xs(t) % value(XS_TOTAL,k) = this % xs(t) % &
value(XS_TOTAL,k) + rx % xs(t, k - j + 1)
end do
! Calculate photon production cross section
do k = 1, size(rx % products)
if (rx % products(k) % particle == PHOTON) then
do k = 1, rx % products_size()
if (rx % product_particle(k) == PHOTON) then
do l = 1, n
this % xs(t) % value(XS_PHOTON_PROD,l+j-1) = &
this % xs(t) % value(XS_PHOTON_PROD,l+j-1) + &
rx % xs(t) % value(l) * rx % products(k) % &
yield % evaluate(this % grid(t) % energy(l+j-1))
rx % xs(t, l) * rx % product_yield(k, &
this % grid(t) % energy(l+j-1))
end do
end if
end do
! Add contribution to absorption cross section
if (is_disappearance(rx % MT)) then
this % xs(t) % value(XS_ABSORPTION,j:j+n-1) = this % xs(t) % &
value(XS_ABSORPTION,j:j+n-1) + rx % xs(t) % value
do k = j, j + n - 1
this % xs(t) % value(XS_ABSORPTION,k) = this % xs(t) % &
value(XS_ABSORPTION,k) + rx % xs(t, k - j + 1)
end do
end if
! Information about fission reactions
@ -680,39 +684,20 @@ contains
! Add contribution to fission cross section
if (is_fission(rx % MT)) then
this % fissionable = .true.
this % xs(t) % value(XS_FISSION,j:j+n-1) = this % xs(t) % &
value(XS_FISSION,j:j+n-1) + rx % xs(t) % value
do k = j, j + n - 1
this % xs(t) % value(XS_FISSION,k) = this % xs(t) % &
value(XS_FISSION,k) + rx % xs(t, k - j + 1)
! Also need to add fission cross sections to absorption
this % xs(t) % value(XS_ABSORPTION,j:j+n-1) = this % xs(t) % &
value(XS_ABSORPTION,j:j+n-1) + rx % xs(t) % value
! Also need to add fission cross sections to absorption
this % xs(t) % value(XS_ABSORPTION,k) = this % xs(t) % &
value(XS_ABSORPTION,k) + rx % xs(t, k - j + 1)
end do
! Keep track of this reaction for easy searching later
if (t == 1) then
i_fission = i_fission + 1
this % index_fission(i_fission) = i
this % n_fission = this % n_fission + 1
! <<<<<<<<<<<<<<<<<<<<<<<<<<<< REMOVE THIS <<<<<<<<<<<<<<<<<<<<<<<<<
! Before the secondary distribution refactor, when the angle/energy
! distribution was uncorrelated, no angle was actually sampled. With
! the refactor, an angle is always sampled for an uncorrelated
! distribution even when no angle distribution exists in the ACE file
! (isotropic is assumed). To preserve the RNG stream, we explicitly
! mark fission reactions so that we avoid the angle sampling.
do k = 1, size(rx % products)
if (rx % products(k) % particle == NEUTRON) then
do m = 1, size(rx % products(k) % distribution)
associate (aedist => rx % products(k) % distribution(m) % obj)
select type (aedist)
type is (UncorrelatedAngleEnergy)
aedist % fission = .true.
end select
end associate
end do
end if
end do
! <<<<<<<<<<<<<<<<<<<<<<<<<<<< REMOVE THIS <<<<<<<<<<<<<<<<<<<<<<<<<
end if
end if ! fission
end do ! temperature
@ -721,12 +706,13 @@ contains
! Determine number of delayed neutron precursors
if (this % fissionable) then
do i = 1, size(this % reactions(this % index_fission(1)) % products)
if (this % reactions(this % index_fission(1)) % products(i) % &
emission_mode == EMISSION_DELAYED) then
this % n_precursor = this % n_precursor + 1
end if
end do
associate (rx => this % reactions(this % index_fission(1)))
do i = 1, rx % products_size()
if (rx % product_emission_mode(i) == EMISSION_DELAYED) then
this % n_precursor = this % n_precursor + 1
end if
end do
end associate
end if
! Calculate nu-fission cross section
@ -761,36 +747,30 @@ contains
select case (emission_mode)
case (EMISSION_PROMPT)
associate (product => this % reactions(this % index_fission(1)) % products(1))
nu = product % yield % evaluate(E)
associate (rx => this % reactions(this % index_fission(1)))
nu = rx % product_yield(1, E)
end associate
case (EMISSION_DELAYED)
if (this % n_precursor > 0) then
if (present(group) .and. group < &
size(this % reactions(this % index_fission(1)) % products)) then
! If delayed group specified, determine yield immediately
associate(p => this % reactions(this % index_fission(1)) % products(1 + group))
nu = p % yield % evaluate(E)
end associate
associate(rx => this % reactions(this % index_fission(1)))
if (present(group) .and. group < rx % products_size()) then
! If delayed group specified, determine yield immediately
nu = rx % product_yield(1 + group, E)
else
nu = ZERO
else
nu = ZERO
do i = 2, rx % products_size()
! Skip any non-neutron products
if (rx % product_particle(i) /= NEUTRON) exit
associate (rx => this % reactions(this % index_fission(1)))
do i = 2, size(rx % products)
associate (product => rx % products(i))
! Skip any non-neutron products
if (product % particle /= NEUTRON) exit
! Evaluate yield
if (product % emission_mode == EMISSION_DELAYED) then
nu = nu + product % yield % evaluate(E)
end if
end associate
! Evaluate yield
if (rx % product_emission_mode(i) == EMISSION_DELAYED) then
nu = nu + rx % product_yield(i, E)
end if
end do
end associate
end if
end if
end associate
else
nu = ZERO
end if
@ -799,8 +779,8 @@ contains
if (allocated(this % total_nu)) then
nu = this % total_nu % evaluate(E)
else
associate (product => this % reactions(this % index_fission(1)) % products(1))
nu = product % yield % evaluate(E)
associate (rx => this % reactions(this % index_fission(1)))
nu = rx % product_yield(1, E)
end associate
end if
end select
@ -868,6 +848,7 @@ contains
integer :: i_high ! upper logarithmic mapping index
integer :: i_rxn ! reaction index
integer :: j ! index in DEPLETION_RX
integer :: threshold ! threshold energy index
real(8) :: f ! interp factor on nuclide energy grid
real(8) :: kT ! temperature in eV
real(8) :: sig_t, sig_a, sig_f ! Intermediate multipole variables
@ -1009,10 +990,11 @@ contains
! need to specifically check its threshold index
i_rxn = this % reaction_index(DEPLETION_RX(1))
if (i_rxn > 0) then
associate (xs => this % reactions(i_rxn) % xs(i_temp))
associate (rx => this % reactions(i_rxn))
threshold = rx % xs_threshold(i_temp)
micro_xs % reaction(1) = (ONE - f) * &
xs % value(i_grid - xs % threshold + 1) + &
f * xs % value(i_grid - xs % threshold + 2)
rx % xs(i_temp, i_grid - threshold + 1) + &
f * rx % xs(i_temp, i_grid - threshold + 2)
end associate
end if
@ -1022,11 +1004,12 @@ contains
! reaction xs appropriately
i_rxn = this % reaction_index(DEPLETION_RX(j))
if (i_rxn > 0) then
associate (xs => this % reactions(i_rxn) % xs(i_temp))
if (i_grid >= xs % threshold) then
associate (rx => this % reactions(i_rxn))
threshold = rx % xs_threshold(i_temp)
if (i_grid >= threshold) then
micro_xs % reaction(j) = (ONE - f) * &
xs % value(i_grid - xs % threshold + 1) + &
f * xs % value(i_grid - xs % threshold + 2)
rx % xs(i_temp, i_grid - threshold + 1) + &
f * rx % xs(i_temp, i_grid - threshold + 2)
elseif (j >= 4) then
! One can show that the the threshold for (n,(x+1)n) is always
! higher than the threshold for (n,xn). Thus, if we are below
@ -1091,8 +1074,9 @@ contains
f = micro_xs % interp_factor
if (i_temp > 0) then
associate (xs => this % reactions(1) % xs(i_temp) % value)
micro_xs % elastic = (ONE - f) * xs(i_grid) + f * xs(i_grid + 1)
associate (rx => this % reactions(1))
micro_xs % elastic = (ONE - f) * rx % xs(i_temp, i_grid) + &
f * rx % xs(i_temp, i_grid + 1)
end associate
else
! For multipole, elastic is total - absorption
@ -1459,6 +1443,7 @@ contains
integer :: i_energy ! index for energy
integer :: i_low ! band index at lower bounding energy
integer :: i_up ! band index at upper bounding energy
integer :: threshold ! threshold energy index
real(8) :: f ! interpolation factor
real(8) :: r ! pseudo-random number
real(8) :: elastic ! elastic cross section
@ -1551,10 +1536,11 @@ contains
f = micro_xs % interp_factor
! Determine inelastic scattering cross section
associate (xs => this % reactions(this % urr_inelastic) % xs(i_temp))
if (i_energy >= xs % threshold) then
inelastic = (ONE - f) * xs % value(i_energy - xs % threshold + 1) + &
f * xs % value(i_energy - xs % threshold + 2)
associate (rx => this % reactions(this % urr_inelastic))
threshold = rx % xs_threshold(i_temp)
if (i_energy >= threshold) then
inelastic = (ONE - f) * rx % xs(i_temp, i_energy - threshold + 1) + &
f * rx % xs(i_temp, i_energy - threshold + 2)
end if
end associate
end if

56
src/physics.F90
View file

@ -506,6 +506,7 @@ contains
integer :: i
integer :: i_grid
integer :: i_temp
integer :: threshold
real(8) :: f
real(8) :: prob
real(8) :: cutoff
@ -545,13 +546,14 @@ contains
FISSION_REACTION_LOOP: do i = 1, nuc % n_fission
i_reaction = nuc % index_fission(i)
associate (xs => nuc % reactions(i_reaction) % xs(i_temp))
associate (rx => nuc % reactions(i_reaction))
! if energy is below threshold for this reaction, skip it
if (i_grid < xs % threshold) cycle
threshold = rx % xs_threshold(i_temp)
if (i_grid < threshold) cycle
! add to cumulative probability
prob = prob + ((ONE - f) * xs % value(i_grid - xs % threshold + 1) &
+ f*(xs % value(i_grid - xs % threshold + 2)))
prob = prob + ((ONE - f) * rx % xs(i_temp, i_grid - threshold + 1) &
+ f*(rx % xs(i_temp, i_grid - threshold + 2)))
end associate
! Create fission bank sites if fission occurs
@ -593,17 +595,17 @@ contains
! Loop through each reaction type
REACTION_LOOP: do i_reaction = 1, size(nuc % reactions)
associate (rx => nuc % reactions(i_reaction))
threshold = rx % xs(i_temp) % threshold
threshold = rx % xs_threshold(i_temp)
! if energy is below threshold for this reaction, skip it
if (i_grid < threshold) cycle
do i_product = 1, size(rx % products)
if (rx % products(i_product) % particle == PHOTON) then
do i_product = 1, rx % products_size()
if (rx % product_particle(i_product) == PHOTON) then
! add to cumulative probability
yield = rx % products(i_product) % yield % evaluate(E)
prob = prob + ((ONE - f) * rx % xs(i_temp) % value(i_grid - threshold + 1) &
+ f*(rx % xs(i_temp) % value(i_grid - threshold + 2))) * yield
yield = rx % product_yield(i_product, E)
prob = prob + ((ONE - f) * rx % xs(i_temp, i_grid - threshold + 1) &
+ f*(rx % xs(i_temp, i_grid - threshold + 2))) * yield
if (prob > cutoff) return
last_valid_reaction = i_reaction
@ -672,6 +674,7 @@ contains
integer :: j
integer :: i_temp
integer :: i_grid
integer :: threshold
real(8) :: f
real(8) :: prob
real(8) :: cutoff
@ -750,14 +753,14 @@ contains
&// trim(nuc % name))
end if
associate (rx => nuc % reactions(i), &
xs => nuc % reactions(i) % xs(i_temp))
associate (rx => nuc % reactions(i))
! if energy is below threshold for this reaction, skip it
if (i_grid < xs % threshold) cycle
threshold = rx % xs_threshold(i_temp)
if (i_grid < threshold) cycle
! add to cumulative probability
prob = prob + ((ONE - f)*xs % value(i_grid - xs % threshold + 1) &
+ f*(xs % value(i_grid - xs % threshold + 2)))
prob = prob + ((ONE - f)*rx % xs(i_temp, i_grid - threshold + 1) &
+ f*(rx % xs(i_temp, i_grid - threshold + 2)))
end associate
end do
@ -839,14 +842,7 @@ contains
vel = sqrt(dot_product(v_n, v_n))
! Sample scattering angle
select type (dist => rxn % products(1) % distribution(1) % obj)
type is (UncorrelatedAngleEnergy)
if (allocated(dist % angle % energy)) then
mu_cm = dist % angle % sample(E)
else
mu_cm = TWO*prn() - ONE
end if
end select
mu_cm = rxn % sample_elastic_mu(E)
! Determine direction cosines in CM
uvw_cm = v_n/vel
@ -1581,7 +1577,7 @@ contains
do group = 1, nuc % n_precursor
! determine delayed neutron precursor yield for group j
yield = rxn % products(1 + group) % yield % evaluate(E_in)
yield = rxn % product_yield(1 + group, E_in)
! Check if this group is sampled
prob = prob + yield
@ -1600,7 +1596,7 @@ contains
do
! sample from energy/angle distribution -- note that mu has already been
! sampled above and doesn't need to be resampled
call rxn % products(1 + group) % sample(E_in, site % E, mu)
call rxn % product_sample(1 + group, E_in, site % E, mu)
! resample if energy is greater than maximum neutron energy
if (site % E < energy_max(NEUTRON)) exit
@ -1624,7 +1620,7 @@ contains
! sample from prompt neutron energy distribution
n_sample = 0
do
call rxn % products(1) % sample(E_in, site % E, mu)
call rxn % product_sample(1, E_in, site % E, mu)
! resample if energy is greater than maximum neutron energy
if (site % E < energy_max(NEUTRON)) exit
@ -1663,7 +1659,7 @@ contains
E_in = p % E
! sample outgoing energy and scattering cosine
call rxn % products(1) % sample(E_in, E, mu)
call rxn % product_sample(1, E_in, E, mu)
! if scattering system is in center-of-mass, transfer cosine of scattering
! angle and outgoing energy from CM to LAB
@ -1692,7 +1688,7 @@ contains
p % coord(1) % uvw = rotate_angle(p % coord(1) % uvw, mu)
! evaluate yield
yield = rxn % products(1) % yield % evaluate(E_in)
yield = rxn % product_yield(1, E_in)
if (mod(yield, ONE) == ZERO) then
! If yield is integral, create exactly that many secondary particles
do i = 1, nint(yield) - 1
@ -1740,8 +1736,8 @@ contains
call sample_photon_product(i_nuclide, p % E, i_reaction, i_product)
! Sample the outgoing energy and angle
call nuclides(i_nuclide) % reactions(i_reaction) % products(i_product) &
% sample(p % E, E, mu)
call nuclides(i_nuclide) % reactions(i_reaction) % &
product_sample(i_product, p % E, E, mu)
! Sample the new direction
uvw = rotate_angle(p % coord(1) % uvw, mu)

115
src/reaction.cpp
View file

@ -4,6 +4,9 @@
#include <utility> // for move
#include "hdf5_interface.h"
#include "endf.h"
#include "random_lcg.h"
#include "secondary_uncorrelated.h"
namespace openmc {
@ -43,6 +46,116 @@ Reaction::Reaction(hid_t group, const std::vector<int>& temperatures)
close_group(pgroup);
}
}
// <<<<<<<<<<<<<<<<<<<<<<<<<<<< REMOVE THIS <<<<<<<<<<<<<<<<<<<<<<<<<
// Before the secondary distribution refactor, when the angle/energy
// distribution was uncorrelated, no angle was actually sampled. With
// the refactor, an angle is always sampled for an uncorrelated
// distribution even when no angle distribution exists in the ACE file
// (isotropic is assumed). To preserve the RNG stream, we explicitly
// mark fission reactions so that we avoid the angle sampling.
if (is_fission(mt_)) {
for (auto& p : products_) {
if (p.particle_ == ParticleType::neutron) {
for (auto& d : p.distribution_) {
auto d_ = dynamic_cast<UncorrelatedAngleEnergy*>(d.get());
if (d_) d_->fission_ = true;
}
}
}
}
// <<<<<<<<<<<<<<<<<<<<<<<<<<<< REMOVE THIS <<<<<<<<<<<<<<<<<<<<<<<<<
}
}
//==============================================================================
// Fortran compatibility functions
//==============================================================================
Reaction* reaction_from_hdf5(hid_t group, int* temperatures, int n)
{
std::vector<int> temps {temperatures, temperatures + n};
return new Reaction{group, temps};
}
void reaction_delete(Reaction* rx) { delete rx; }
int reaction_mt(Reaction* rx) { return rx->mt_; }
double reaction_q_value(Reaction* rx) { return rx->q_value_; }
bool reaction_scatter_in_cm(Reaction* rx) { return rx->scatter_in_cm_; }
double reaction_product_decay_rate(Reaction* rx, int product)
{
return rx->products_[product - 1].decay_rate_;
}
int reaction_product_emission_mode(Reaction* rx, int product)
{
switch (rx->products_[product - 1].emission_mode_) {
case ReactionProduct::EmissionMode::prompt:
return 1;
case ReactionProduct::EmissionMode::delayed:
return 2;
case ReactionProduct::EmissionMode::total:
return 3;
}
}
int reaction_product_particle(Reaction* rx, int product)
{
switch (rx->products_[product - 1].particle_) {
case ParticleType::neutron:
return 1;
case ParticleType::photon:
return 2;
case ParticleType::electron:
return 3;
case ParticleType::positron:
return 4;
}
}
void reaction_product_sample(Reaction* rx, int product, double E_in, double* E_out, double* mu)
{
rx->products_[product - 1].sample(E_in, *E_out, *mu);
}
double reaction_product_yield(Reaction* rx, int product, double E)
{
return (*rx->products_[product - 1].yield_)(E);
}
int reaction_products_size(Reaction* rx) { return rx->products_.size(); }
double reaction_xs(Reaction* rx, int temperature, int energy)
{
return rx->xs_[temperature - 1].value[energy - 1];
}
double reaction_sample_elastic_mu(Reaction* rx, double E)
{
// Get elastic scattering distribution
auto& d = rx->products_[0].distribution_[0];
// Check if it is an uncorrelated angle-energy distribution
auto d_ = dynamic_cast<UncorrelatedAngleEnergy*>(d.get());
if (d_) {
return d_->sampleMu(E);
} else {
return 2.0*prn() - 1.0;
}
}
int reaction_xs_size(Reaction* rx, int temperature)
{
return rx->xs_[temperature - 1].value.size();
}
int reaction_xs_threshold(Reaction* rx, int temperature)
{
return rx->xs_[temperature - 1].threshold;
}
}

23
src/reaction.h
View file

@ -24,6 +24,29 @@ public:
std::vector<ReactionProduct> products_;
};
//==============================================================================
// Fortran compatibility functions
//==============================================================================
extern "C" {
Reaction* reaction_from_hdf5(hid_t group, int* temperatures, int n);
void reaction_delete(Reaction* rx);
int reaction_mt(Reaction* rx);
double reaction_q_value(Reaction* rx);
bool reaction_scatter_in_cm(Reaction* rx);
double reaction_product_decay_rate(Reaction* rx, int product);
int reaction_product_emission_mode(Reaction* rx, int product);
int reaction_product_particle(Reaction* rx, int product);
void reaction_product_sample(Reaction* rx, int product, double E_in,
double* E_out, double* mu);
int reaction_products_size(Reaction* rx);
double reaction_product_yield(Reaction* rx, int product, double E);
double reaction_sample_elastic_mu(Reaction* rx, double E);
double reaction_xs(Reaction* xs, int temperature, int energy);
int reaction_xs_size(Reaction* xs, int temperature);
int reaction_xs_threshold(Reaction* xs, int temperature);
}
}
#endif // OPENMC_REACTION_H

293
src/reaction_header.F90
View file

@ -1,83 +1,268 @@
module reaction_header
use, intrinsic :: ISO_C_BINDING
use constants, only: MAX_WORD_LEN
use hdf5_interface
use product_header, only: ReactionProduct
use stl_vector, only: VectorInt
use string, only: to_str, starts_with
implicit none
private
!===============================================================================
! REACTION contains the cross-section and secondary energy and angle
! distributions for a single reaction in a continuous-energy ACE-format table
!===============================================================================
type TemperatureXS
integer :: threshold ! Energy grid index of threshold
real(8), allocatable :: value(:) ! Cross section values
end type TemperatureXS
type Reaction
integer :: MT ! ENDF MT value
real(8) :: Q_value ! Reaction Q value
logical :: scatter_in_cm ! scattering system in center-of-mass?
type(TemperatureXS), allocatable :: xs(:)
type(ReactionProduct), allocatable :: products(:)
type, public :: Reaction
type(C_PTR) :: ptr
integer(C_INT) :: MT ! ENDF MT value
real(C_DOUBLE) :: Q_value ! Reaction Q value
logical(C_BOOL) :: scatter_in_cm ! scattering system in center-of-mass?
contains
procedure :: from_hdf5 => reaction_from_hdf5
procedure :: from_hdf5
procedure :: mt_
procedure :: q_value_
procedure :: scatter_in_cm_
procedure :: product_decay_rate
procedure :: product_emission_mode
procedure :: product_particle
procedure :: product_sample
procedure :: product_yield
procedure :: products_size
procedure :: sample_elastic_mu
procedure :: xs
procedure :: xs_size
procedure :: xs_threshold
end type Reaction
interface
function reaction_from_hdf5(group, temperatures, n) result(ptr) bind(C)
import C_PTR, HID_T, C_INT
integer(HID_T), value :: group
integer(C_INT), intent(in) :: temperatures
integer(C_INT), value :: n
type(C_PTR) :: ptr
end function
function reaction_mt(ptr) result(mt) bind(C)
import C_PTR, C_INT
type(C_PTR), value :: ptr
integer(C_INT) :: mt
end function
function reaction_q_value(ptr) result(q_value) bind(C)
import C_PTR, C_DOUBLE
type(C_PTR), value :: ptr
real(C_DOUBLE) :: q_value
end function
function reaction_scatter_in_cm(ptr) result(b) bind(C)
import C_PTR, C_BOOL
type(C_PTR), value :: ptr
logical(C_BOOL) :: b
end function
pure function reaction_product_decay_rate(ptr, product) result(rate) bind(C)
import C_PTR, C_INT, C_DOUBLE
type(C_PTR), value :: ptr
integer(C_INT), value :: product
real(C_DOUBLE) :: rate
end function
pure function reaction_product_emission_mode(ptr, product) result(m) bind(C)
import C_PTR, C_INT
type(C_PTR), value :: ptr
integer(C_INT), value :: product
integer(C_INT) :: m
end function
pure function reaction_product_particle(ptr, product) result(particle) bind(C)
import C_PTR, C_INT
type(C_PTR), value :: ptr
integer(C_INT), value :: product
integer(C_INT) :: particle
end function
subroutine reaction_product_sample(ptr, product, E_in, E_out, mu) bind(C)
import C_PTR, C_INT, C_DOUBLE
type(C_PTR), value :: ptr
integer(C_INT), value :: product
real(C_DOUBLE), value :: E_in
real(C_DOUBLE), intent(out) :: E_out
real(C_DOUBLE), intent(out) :: mu
end subroutine
pure function reaction_product_yield(ptr, product, E) result(val) bind(C)
import C_PTR, C_INT, C_DOUBLE
type(C_PTR), value :: ptr
integer(C_INT), value :: product
real(C_DOUBLE), value :: E
real(C_DOUBLE) :: val
end function
pure function reaction_products_size(ptr) result(sz) bind(C)
import C_PTR, C_INT
type(C_PTR), value :: ptr
integer(C_INT) :: sz
end function
function reaction_sample_elastic_mu(ptr, E) result(mu) bind(C)
import C_PTR, C_INT, C_DOUBLE
type(C_PTR), value :: ptr
real(C_DOUBLE), value :: E
real(C_DOUBLE) :: mu
end function
function reaction_xs(ptr, temperature, energy) result(xs) bind(C)
import C_PTR, C_INT, C_DOUBLE
type(C_PTR), value :: ptr
integer(C_INT), value :: temperature
integer(C_INT), value :: energy
real(C_DOUBLE) :: xs
end function
function reaction_xs_size(ptr, temperature) result(sz) bind(C)
import C_PTR, C_INT
type(C_PTR), value :: ptr
integer(C_INT), value :: temperature
integer(C_INT) :: sz
end function
function reaction_xs_threshold(ptr, temperature) result(threshold) bind(C)
import C_PTR, C_INT
type(C_PTR), value :: ptr
integer(C_INT), value :: temperature
integer(C_INT) :: threshold
end function
end interface
contains
subroutine reaction_from_hdf5(this, group_id, temperatures)
subroutine from_hdf5(this, group_id, temperatures)
class(Reaction), intent(inout) :: this
integer(HID_T), intent(in) :: group_id
type(VectorInt), intent(in) :: temperatures
integer :: i
integer :: cm
integer :: n_product
integer(HID_T) :: pgroup
integer(HID_T) :: xs, temp_group
integer(HSIZE_T) :: dims(1)
integer(HSIZE_T) :: j
character(MAX_WORD_LEN) :: temp_str ! temperature dataset name, e.g. '294K'
character(MAX_WORD_LEN), allocatable :: grp_names(:)
integer(C_INT) :: dummy
integer(C_INT) :: n
call read_attribute(this % Q_value, group_id, 'Q_value')
call read_attribute(this % MT, group_id, 'mt')
call read_attribute(cm, group_id, 'center_of_mass')
this % scatter_in_cm = (cm == 1)
n = temperatures % size()
if (n > 0) then
this % ptr = reaction_from_hdf5(group_id, temperatures % data(1), n)
else
! In this case, temperatures % data(1) doesn't exist, so we just pass a
! dummy value
this % ptr = reaction_from_hdf5(group_id, dummy, n)
end if
this % MT = reaction_mt(this % ptr)
this % Q_value = reaction_q_value(this % ptr)
this % scatter_in_cm = reaction_scatter_in_cm(this % ptr)
end subroutine from_hdf5
! Read cross section and threshold_idx data
allocate(this % xs(temperatures % size()))
do i = 1, temperatures % size()
temp_str = trim(to_str(temperatures % data(i))) // "K"
temp_group = open_group(group_id, temp_str)
xs = open_dataset(temp_group, 'xs')
call read_attribute(this % xs(i) % threshold, xs, 'threshold_idx')
call get_shape(xs, dims)
allocate(this % xs(i) % value(dims(1)))
call read_dataset(this % xs(i) % value, xs)
call close_dataset(xs)
call close_group(temp_group)
end do
function mt_(this) result(mt)
class(Reaction), intent(in) :: this
integer(C_INT) :: MT
! Determine number of products
n_product = 0
call get_groups(group_id, grp_names)
do j = 1, size(grp_names)
if (starts_with(grp_names(j), "product_")) n_product = n_product + 1
end do
mt = reaction_mt(this % ptr)
end function
! Read products
allocate(this % products(n_product))
do i = 1, n_product
pgroup = open_group(group_id, 'product_' // trim(to_str(i - 1)))
call this % products(i) % from_hdf5(pgroup)
call close_group(pgroup)
end do
end subroutine reaction_from_hdf5
function q_value_(this) result(q_value)
class(Reaction), intent(in) :: this
real(C_DOUBLE) :: q_value
q_value = reaction_q_value(this % ptr)
end function
function scatter_in_cm_(this) result(cm)
class (Reaction), intent(in) :: this
logical(C_BOOL) :: cm
cm = reaction_scatter_in_cm(this % ptr)
end function
pure function product_decay_rate(this, product) result(rate)
class(Reaction), intent(in) :: this
integer(C_INT), intent(in) :: product
real(C_DOUBLE) :: rate
rate = reaction_product_decay_rate(this % ptr, product)
end function
pure function product_emission_mode(this, product) result(m)
class(Reaction), intent(in) :: this
integer(C_INT), intent(in) :: product
integer(C_INT) :: m
m = reaction_product_emission_mode(this % ptr, product)
end function
pure function product_particle(this, product) result(p)
class(Reaction), intent(in) :: this
integer(C_INT), intent(in) :: product
integer(C_INT) :: p
p = reaction_product_particle(this % ptr, product)
end function
subroutine product_sample(this, product, E_in, E_out, mu)
class(Reaction), intent(in) :: this
integer(C_INT), intent(in) :: product
real(C_DOUBLE), intent(in) :: E_in
real(C_DOUBLE), intent(out) :: E_out
real(C_DOUBLE), intent(out) :: mu
call reaction_product_sample(this % ptr, product, E_in, E_out, mu)
end subroutine
pure function product_yield(this, product, E) result(val)
class(Reaction), intent(in) :: this
integer(C_INT), intent(in) :: product
real(C_DOUBLE), intent(in) :: E
real(C_DOUBLE) :: val
val = reaction_product_yield(this % ptr, product, E)
end function
pure function products_size(this) result(sz)
class(Reaction), intent(in) :: this
integer(C_INT) :: sz
sz = reaction_products_size(this % ptr)
end function
function sample_elastic_mu(this, E) result(mu)
class(Reaction), intent(in) :: this
real(C_DOUBLE), intent(in) :: E
real(C_DOUBLE) :: mu
mu = reaction_sample_elastic_mu(this % ptr, E)
end function
function xs(this, temperature, energy) result(val)
class(Reaction), intent(in) :: this
integer(C_INT), intent(in) :: temperature
integer(C_INT), intent(in) :: energy
real(C_DOUBLE) :: val
val = reaction_xs(this % ptr, temperature, energy)
end function
function xs_size(this, temperature) result(sz)
class(Reaction), intent(in) :: this
integer(C_INT) :: temperature
integer(C_INT) :: sz
sz = reaction_xs_size(this % ptr, temperature)
end function
function xs_threshold(this, temperature) result(val)
class(Reaction), intent(in) :: this
integer(C_INT), intent(in) :: temperature
integer(C_INT) :: val
val = reaction_xs_threshold(this % ptr, temperature)
end function
end module reaction_header

2
src/reaction_product.h
View file

@ -15,7 +15,7 @@ class ReactionProduct {
public:
explicit ReactionProduct(hid_t group);
void sample(double E_in, double& E_out, double& mu) const;
private:
enum class EmissionMode {
prompt, // Prompt emission of secondary particle
total, // Delayed emission of secondary particle

4
src/secondary_uncorrelated.h
View file

@ -15,8 +15,10 @@ class UncorrelatedAngleEnergy : public AngleEnergy {
public:
explicit UncorrelatedAngleEnergy(hid_t group);
void sample(double E_in, double& E_out, double& mu) const;
private:
double sampleMu(double E) const { return angle_.sample(E); }
bool fission_ {false};
private:
AngleDistribution angle_;
std::unique_ptr<EnergyDistribution> energy_;
};

49
src/tallies/tally.F90
View file

@ -94,6 +94,7 @@ contains
integer :: k ! loop index for bank sites
integer :: d_bin ! delayed group bin index
integer :: dg_filter ! index of delayed group filter
integer :: threshold ! threshold energy index
real(8) :: yield ! delayed neutron yield
real(8) :: atom_density_ ! atom/b-cm
real(8) :: f ! interpolation factor
@ -227,7 +228,7 @@ contains
! Get yield and apply to score
associate (rxn => nuclides(p % event_nuclide) % reactions(m))
score = p % last_wgt * flux &
* rxn % products(1) % yield % evaluate(E)
* rxn % product_yield(1, E)
end associate
end if
@ -633,7 +634,7 @@ contains
score = p % absorb_wgt * yield * &
micro_xs(p % event_nuclide) % fission &
/ micro_xs(p % event_nuclide) % absorption &
* rxn % products(1 + d) % decay_rate * flux
* rxn % product_decay_rate(1 + d) * flux
end associate
! Tally to bin
@ -657,9 +658,8 @@ contains
! rxn % products array to be exceeded. Hence, we use the size
! of this array and not the MAX_DELAYED_GROUPS constant for
! this loop.
do d = 1, size(rxn % products) - 2
score = score + rxn % products(1 + d) % decay_rate * &
do d = 1, rxn % products_size() - 2
score = score + rxn % product_decay_rate(1 + d) * &
p % absorb_wgt &
* micro_xs(p % event_nuclide) % fission &
* nuclides(p % event_nuclide) % &
@ -699,7 +699,7 @@ contains
! determine score based on bank site weight and keff.
score = score + keff * fission_bank(n_bank - p % n_bank + k) &
% wgt * rxn % products(1 + g) % decay_rate * flux
% wgt * rxn % product_decay_rate(1 + g) * flux
end associate
! if the delayed group filter is present, tally to corresponding
@ -755,7 +755,7 @@ contains
! Compute the score and tally to bin
score = micro_xs(i_nuclide) % fission * yield * flux * &
atom_density * rxn % products(1 + d) % decay_rate
atom_density * rxn % product_decay_rate(1 + d)
end associate
! Tally to bin
@ -778,11 +778,10 @@ contains
! groups since this could cause the range of the rxn % products
! array to be exceeded. Hence, we use the size of this array
! and not the MAX_DELAYED_GROUPS constant for this loop.
do d = 1, size(rxn % products) - 2
do d = 1, rxn % products_size() - 2
score = score + micro_xs(i_nuclide) % fission * flux * &
nuclides(i_nuclide) % nu(E, EMISSION_DELAYED) * &
atom_density * rxn % products(1 + d) % decay_rate
atom_density * rxn % product_decay_rate(1 + d)
end do
end associate
end if
@ -824,7 +823,7 @@ contains
! Compute the score
score = micro_xs(i_nuc) % fission * yield * flux * &
atom_density_ &
* rxn % products(1 + d) % decay_rate
* rxn % product_decay_rate(1 + d)
end associate
! Tally to bin
@ -860,13 +859,13 @@ contains
! rxn % products array to be exceeded. Hence, we use the
! size of this array and not the MAX_DELAYED_GROUPS
! constant for this loop.
do d = 1, size(rxn % products) - 2
do d = 1, rxn % products_size() - 2
! Accumulate the contribution from each nuclide
score = score + micro_xs(i_nuc) % fission &
* nuclides(i_nuc) % nu(E, EMISSION_DELAYED) &
* atom_density_ * flux &
* rxn % products(1 + d) % decay_rate
* rxn % product_decay_rate(1 + d)
end do
end associate
end if
@ -1132,12 +1131,12 @@ contains
i_energy = micro_xs(i_nuclide) % index_grid
f = micro_xs(i_nuclide) % interp_factor
associate (xs => nuclides(i_nuclide) % reactions(m) &
% xs(i_temp))
if (i_energy >= xs % threshold) then
score = ((ONE - f) * xs % value(i_energy - &
xs % threshold + 1) + f * xs % value(i_energy - &
xs % threshold + 2)) * atom_density * flux
associate (rx => nuclides(i_nuclide) % reactions(m))
threshold = rx % xs_threshold(i_temp)
if (i_energy >= threshold) then
score = ((ONE - f) * rx % xs(i_temp, i_energy - &
threshold + 1) + f * rx % xs(i_temp, i_energy - &
threshold + 2)) * atom_density * flux
end if
end associate
else
@ -1165,12 +1164,12 @@ contains
i_energy = micro_xs(i_nuc) % index_grid
f = micro_xs(i_nuc) % interp_factor
associate (xs => nuclides(i_nuc) % reactions(m) &
% xs(i_temp))
if (i_energy >= xs % threshold) then
score = score + ((ONE - f) * xs % value(i_energy - &
xs % threshold + 1) + f * xs % value(i_energy - &
xs % threshold + 2)) * atom_density_ * flux
associate (rx => nuclides(i_nuc) % reactions(m))
threshold = rx % xs_threshold(i_temp)
if (i_energy >= threshold) then
score = score + ((ONE - f) * rx % xs(i_temp, i_energy - &
threshold + 1) + f * rx % xs(i_temp, i_energy - &
threshold + 2)) * atom_density_ * flux
end if
end associate
else