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

Convert sample_fission and sample_fission_neutron to C++

Browse source
This commit is contained in:
Paul Romano 2018-11-19 07:47:27 -06:00
parent 934e26ca53
commit 1f90ec1560
5 changed files with 222 additions and 312 deletions
Download patch file Download diff file Expand all files Collapse all files

14
include/openmc/nuclide.h
View file

@ -11,7 +11,9 @@
#include <hdf5.h>
#include "openmc/constants.h"
#include "openmc/endf.h"
#include "openmc/reaction.h"
#include "openmc/reaction_product.h"
namespace openmc {
@ -21,16 +23,28 @@ namespace openmc {
class Nuclide {
public:
using EmissionMode = ReactionProduct::EmissionMode;
// Constructors
Nuclide(hid_t group, const double* temperature, int n);
// Methods
double nu(double E, EmissionMode mode, int group=0);
// Data members
std::string name_; //! Name of nuclide, e.g. "U235"
int Z_; //! Atomic number
int A_; //! Mass number
int metastable_; //! Metastable state
double awr_; //! Atomic weight ratio
std::vector<double> kTs_; //! temperatures in eV (k*T)
bool fissionable_ {false}; //! Whether nuclide is fissionable
bool has_partial_fission_ {false}; //! has partial fission reactions?
std::vector<Reaction*> fission_rx_; //! Fission reactions
int n_precursor_ {0}; //! Number of delayed neutron precursors
std::unique_ptr<Function1D> total_nu_; //! Total neutron yield
std::vector<std::unique_ptr<Reaction>> reactions_; //! Reactions
private:

6
include/openmc/physics.h
View file

@ -40,12 +40,12 @@ void sample_nuclide(const Particle* p, int mt, int* i_nuclide, int* i_nuc_mat);
//! Determine the average total, prompt, and delayed neutrons produced from
//! fission and creates appropriate bank sites.
void create_fission_sites(Particle* p, int i_nuclide, int i_rx,
void create_fission_sites(Particle* p, int i_nuclide, const Reaction* rx,
Bank* bank_array, int64_t* bank_size, int64_t bank_capacity);
// void sample_element(Particle* p);
extern "C" int sample_fission(int i_nuclide, double E);
Reaction* sample_fission(int i_nuclide, double E);
// void sample_photon_product(int i_nuclide, double E, int* i_rx, int* i_product);
@ -64,7 +64,7 @@ extern "C" void scatter(Particle*, int i_nuclide, int i_nuc_mat);
// void sample_cxs_target_velocity(int i_nuclide, Direction* v_target, double E, Direction u,
// double kT);
extern "C" void sample_fission_neutron(int i_nuclide, int i_rx, double E_in, Bank* site);
void sample_fission_neutron(int i_nuclide, const Reaction* rx, double E_in, Bank* site);
// void inelastic_scatter(int i_nuclide, const Reaction& rx, Particle* p);

95
src/nuclide.cpp
View file

@ -52,7 +52,6 @@ Nuclide::Nuclide(hid_t group, const double* temperature, int n)
temps_available.push_back(T / K_BOLTZMANN);
}
std::sort(temps_available.begin(), temps_available.end());
close_group(kT_group);
// If only one temperature is available, revert to nearest temperature
if (temps_available.size() == 1 && settings::temperature_method == TEMPERATURE_INTERPOLATION) {
@ -144,6 +143,15 @@ Nuclide::Nuclide(hid_t group, const double* temperature, int n)
// Sort temperatures to read
std::sort(temps_to_read.begin(), temps_to_read.end());
// Determine exact kT values
for (const auto& T : temps_to_read) {
std::string dset {std::to_string(T) + "K"};
double kT;
read_dataset(kT_group, dset.c_str(), kT);
kTs_.push_back(kT);
}
close_group(kT_group);
// Read reactions
hid_t rxs_group = open_group(group, "reactions");
for (auto name : group_names(rxs_group)) {
@ -155,17 +163,92 @@ Nuclide::Nuclide(hid_t group, const double* temperature, int n)
}
close_group(rxs_group);
// Check for nu-total
if (object_exists(group, "total_nu")) {
// Read total nu data
hid_t nu_group = open_group(group, "total_nu");
hid_t nu_dset = open_dataset(nu_group, "yield");
std::string func_type;
read_attribute(nu_dset, "type", func_type);
if (func_type == "Tabulated1D") {
total_nu_ = std::make_unique<Tabulated1D>(nu_dset);
} else if (func_type == "Polynomial") {
total_nu_ = std::make_unique<Polynomial>(nu_dset);
}
close_dataset(nu_dset);
close_group(nu_group);
}
this->create_derived();
}
void Nuclide::create_derived()
{
for (const auto& rx : reactions_) {
// Skip redundant reactions
if (rx->redundant_) continue;
for (int i = 0; i < reactions_.size(); ++i) {
const auto& rx {reactions_[i]};
if (is_fission(rx->mt_)) {
fissionable_ = true;
for (int t = 0; t < kTs_.size(); ++t) {
// Skip redundant reactions
if (rx->redundant_) continue;
if (is_fission(rx->mt_)) {
fissionable_ = true;
// Keep track of fission reactions
if (t == 0) {
fission_rx_.push_back(rx.get());
if (rx->mt_ == N_F) has_partial_fission_ = true;
}
}
}
}
// Determine number of delayed neutron precursors
if (fissionable_) {
for (const auto& product : fission_rx_[0]->products_) {
if (product.emission_mode_ == EmissionMode::delayed) {
++n_precursor_;
}
}
}
}
double Nuclide::nu(double E, EmissionMode mode, int group)
{
if (!fissionable_) return 0.0;
switch (mode) {
case EmissionMode::prompt:
return (*fission_rx_[0]->products_[0].yield_)(E);
case EmissionMode::delayed:
if (n_precursor_ > 0) {
auto rx = fission_rx_[0];
if (group >= 1 && group < rx->products_.size()) {
// If delayed group specified, determine yield immediately
return (*rx->products_[group].yield_)(E);
} else {
double nu {0.0};
for (int i = 1; i < rx->products_.size(); ++i) {
// Skip any non-neutron products
const auto& product = rx->products_[i];
if (product.particle_ != ParticleType::neutron) continue;
// Evaluate yield
if (product.emission_mode_ == EmissionMode::delayed) {
nu += (*product.yield_)(E);
}
}
return nu;
}
} else {
return 0.0;
}
case EmissionMode::total:
if (total_nu_) {
return (*total_nu_)(E);
} else {
return (*fission_rx_[0]->products_[0].yield_)(E);
}
}
}

180
src/physics.F90
View file

@ -254,74 +254,6 @@ contains
end function sample_element
!===============================================================================
! SAMPLE_FISSION
!===============================================================================
function sample_fission(i_nuclide, E) result (i_reaction) bind(C)
integer(C_INT), value :: i_nuclide ! index in nuclides array
real(C_DOUBLE), value :: E ! incident neutron energy
integer(C_INT) :: i_reaction ! index in nuc % reactions array
integer :: i
integer :: i_grid
integer :: i_temp
integer :: threshold
real(8) :: f
real(8) :: prob
real(8) :: cutoff
type(Nuclide), pointer :: nuc
! Get pointer to nuclide
nuc => nuclides(i_nuclide)
! If we're in the URR, by default use the first fission reaction. We also
! default to the first reaction if we know that there are no partial fission
! reactions
if (micro_xs(i_nuclide) % use_ptable .or. &
.not. nuc % has_partial_fission) then
i_reaction = nuc % index_fission(1)
return
end if
! Check to see if we are in a windowed multipole range. WMP only supports
! the first fission reaction.
if (nuc % mp_present) then
if (E >= nuc % multipole % E_min .and. &
E <= nuc % multipole % E_max) then
i_reaction = nuc % index_fission(1)
return
end if
end if
! Get grid index and interpolatoin factor and sample fission cdf
i_temp = micro_xs(i_nuclide) % index_temp
i_grid = micro_xs(i_nuclide) % index_grid
f = micro_xs(i_nuclide) % interp_factor
cutoff = prn() * micro_xs(i_nuclide) % fission
prob = ZERO
! Loop through each partial fission reaction type
FISSION_REACTION_LOOP: do i = 1, nuc % n_fission
i_reaction = nuc % index_fission(i)
associate (rx => nuc % reactions(i_reaction))
! if energy is below threshold for this reaction, skip it
threshold = rx % xs_threshold(i_temp)
if (i_grid < threshold) cycle
! add to cumulative probability
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
if (prob > cutoff) exit FISSION_REACTION_LOOP
end do FISSION_REACTION_LOOP
end function sample_fission
!===============================================================================
! SAMPLE_PHOTON_PRODUCT
!===============================================================================
@ -928,118 +860,6 @@ contains
end subroutine sample_cxs_target_velocity
!===============================================================================
! SAMPLE_FISSION_NEUTRON
!===============================================================================
subroutine sample_fission_neutron(i_nuc, i_rx, E_in, site) bind(C)
integer(C_INT), value :: i_nuc
integeR(C_INT), value :: i_rx
real(C_DOUBLE), value :: E_in
type(Bank), intent(inout) :: site
integer :: group ! index on nu energy grid / precursor group
integer :: n_sample ! number of resamples
real(8) :: nu_t ! total nu
real(8) :: nu_d ! delayed nu
real(8) :: beta ! delayed neutron fraction
real(8) :: xi ! random number
real(8) :: yield ! delayed neutron precursor yield
real(8) :: prob ! cumulative probability
real(8) :: mu ! cosine of scattering angle
real(8) :: phi ! azimuthal angle
associate (nuc => nuclides(i_nuc), rxn => nuclides(i_nuc) % reactions(i_rx))
! Sample cosine of angle -- fission neutrons are always emitted
! isotropically. Sometimes in ACE data, fission reactions actually have
! an angular distribution listed, but for those that do, it's simply just
! a uniform distribution in mu
mu = TWO * prn() - ONE
! Sample azimuthal angle uniformly in [0,2*pi)
phi = TWO*PI*prn()
site % uvw(1) = mu
site % uvw(2) = sqrt(ONE - mu*mu) * cos(phi)
site % uvw(3) = sqrt(ONE - mu*mu) * sin(phi)
! Determine total nu, delayed nu, and delayed neutron fraction
nu_t = nuc % nu(E_in, EMISSION_TOTAL)
nu_d = nuc % nu(E_in, EMISSION_DELAYED)
beta = nu_d / nu_t
if (prn() < beta) then
! ====================================================================
! DELAYED NEUTRON SAMPLED
! sampled delayed precursor group
xi = prn()*nu_d
prob = ZERO
do group = 1, nuc % n_precursor
! determine delayed neutron precursor yield for group j
yield = rxn % product_yield(1 + group, E_in)
! Check if this group is sampled
prob = prob + yield
if (xi < prob) exit
end do
! if the sum of the probabilities is slightly less than one and the
! random number is greater, j will be greater than nuc %
! n_precursor -- check for this condition
group = min(group, nuc % n_precursor)
! set the delayed group for the particle born from fission
site % delayed_group = group
n_sample = 0
do
! sample from energy/angle distribution -- note that mu has already been
! sampled above and doesn't need to be resampled
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
! check for large number of resamples
n_sample = n_sample + 1
if (n_sample == MAX_SAMPLE) then
! call particle_write_restart(p)
call fatal_error("Resampled energy distribution maximum number of " &
// "times for nuclide " // nuc % name)
end if
end do
else
! ====================================================================
! PROMPT NEUTRON SAMPLED
! set the delayed group for the particle born from fission to 0
site % delayed_group = 0
! sample from prompt neutron energy distribution
n_sample = 0
do
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
! check for large number of resamples
n_sample = n_sample + 1
if (n_sample == MAX_SAMPLE) then
! call particle_write_restart(p)
call fatal_error("Resampled energy distribution maximum number of " &
// "times for nuclide " // nuc % name)
end if
end do
end if
end associate
end subroutine sample_fission_neutron
!===============================================================================
! INELASTIC_SCATTER handles all reactions with a single secondary neutron (other
! than fission), i.e. level scattering, (n,np), (n,na), etc.

239
src/physics.cpp
View file

@ -82,13 +82,13 @@ void sample_neutron_reaction(Particle* p)
const auto& nuc {data::nuclides[i_nuclide-1]};
if (nuc->fissionable_) {
int i_rx = sample_fission(i_nuclide, p->E);
Reaction* rx = sample_fission(i_nuclide, p->E);
if (settings::run_mode == RUN_MODE_EIGENVALUE) {
create_fission_sites(p, i_nuclide, i_rx, simulation::fission_bank.data(),
create_fission_sites(p, i_nuclide, rx, simulation::fission_bank.data(),
&simulation::n_bank, simulation::fission_bank.size());
} else if (settings::run_mode == RUN_MODE_FIXEDSOURCE &&
settings::create_fission_neutrons) {
create_fission_sites(p, i_nuclide, i_rx, p->secondary_bank,
create_fission_sites(p, i_nuclide, rx, p->secondary_bank,
&p->n_secondary, MAX_SECONDARY);
}
}
@ -129,7 +129,7 @@ void sample_neutron_reaction(Particle* p)
}
void
create_fission_sites(Particle* p, int i_nuclide, int i_rx, Bank* bank_array,
create_fission_sites(Particle* p, int i_nuclide, const Reaction* rx, Bank* bank_array,
int64_t* size_bank, int64_t bank_capacity)
{
// TODO: Heat generation from fission
@ -185,7 +185,7 @@ create_fission_sites(Particle* p, int i_nuclide, int i_rx, Bank* bank_array,
bank_array[i].wgt = 1. / weight;
// Sample delayed group and angle/energy for fission reaction
sample_fission_neutron(i_nuclide, i_rx, p->E, &bank_array[i]);
sample_fission_neutron(i_nuclide, rx, p->E, &bank_array[i]);
// Set the delayed group on the particle as well
p->delayed_group = bank_array[i].delayed_group;
@ -485,58 +485,48 @@ void sample_nuclide(const Particle* p, int mt, int* i_nuclide, int* i_nuc_mat)
// end associate
// }
// int sample_fission(int i_nuclide, double E)
// {
// // Get pointer to nuclide
// nuc => nuclides(i_nuclide)
Reaction* sample_fission(int i_nuclide, double E)
{
// Get pointer to nuclide
const auto& nuc {data::nuclides[i_nuclide-1]};
// // If we're in the URR, by default use the first fission reaction. We also
// // default to the first reaction if we know that there are no partial fission
// // reactions
// if (micro_xs(i_nuclide) % use_ptable || &
// !nuc % has_partial_fission) {
// i_reaction = nuc % index_fission(1)
// return
// }
// If we're in the URR, by default use the first fission reaction. We also
// default to the first reaction if we know that there are no partial fission
// reactions
if (simulation::micro_xs[i_nuclide-1].use_ptable || !nuc->has_partial_fission_) {
return nuc->fission_rx_[0];
}
// // Check to see if we are in a windowed multipole range. WMP only supports
// // the first fission reaction.
// if (nuc % mp_present) {
// if (E >= nuc % multipole % E_min && &
// E <= nuc % multipole % E_max) {
// i_reaction = nuc % index_fission(1)
// return
// }
// }
// Check to see if we are in a windowed multipole range. WMP only supports
// the first fission reaction.
// if (nuc % mp_present) {
// if (E >= nuc % multipole % E_min && &
// E <= nuc % multipole % E_max) {
// return nuc->fission_rx_[0];
// }
// }
// // Get grid index and interpolatoin factor and sample fission cdf
// i_temp = micro_xs(i_nuclide) % index_temp
// i_grid = micro_xs(i_nuclide) % index_grid
// f = micro_xs(i_nuclide) % interp_factor
// cutoff = prn() * micro_xs(i_nuclide) % fission
// prob = 0.0
// Get grid index and interpolatoin factor and sample fission cdf
int i_temp = simulation::micro_xs[i_nuclide-1].index_temp;
int i_grid = simulation::micro_xs[i_nuclide-1].index_grid;
double f = simulation::micro_xs[i_nuclide-1].interp_factor;
double cutoff = prn() * simulation::micro_xs[i_nuclide-1].fission;
double prob = 0.0;
// // Loop through each partial fission reaction type
// Loop through each partial fission reaction type
for (auto& rx : nuc->reactions_) {
// if energy is below threshold for this reaction, skip it
int threshold = rx->xs_[i_temp-1].threshold;
if (i_grid < threshold) continue;
// FISSION_REACTION_LOOP: do i = 1, nuc % n_fission
// i_reaction = nuc % index_fission(i)
// add to cumulative probability
prob += (1.0 - f) * rx->xs_[i_temp-1].value[i_grid - threshold]
+ f*rx->xs_[i_temp].value[i_grid - threshold + 1];
// associate (rx => nuc % reactions(i_reaction))
// // if energy is below threshold for this reaction, skip it
// threshold = rx % xs_threshold(i_temp)
// if (i_grid < threshold) cycle
// // add to cumulative probability
// prob = prob + ((1.0 - 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
// if (prob > cutoff) exit FISSION_REACTION_LOOP
// end do FISSION_REACTION_LOOP
// end subroutine sample_fission
// }
// Create fission bank sites if fission occurs
if (prob > cutoff) break;
}
}
// void sample_photon_product(int i_nuclide, double E, int* i_rx, int* i_product);
// {
@ -997,93 +987,96 @@ void sample_nuclide(const Particle* p, int mt, int* i_nuclide, int* i_nuc_mat)
// v_target = vt * rotate_angle(uvw, mu)
// }
// void sample_fission_neutron(int i_nuclide, const Reaction& rx, double E_in, Bank* site)
// {
// // Sample cosine of angle -- fission neutrons are always emitted
// // isotropically. Sometimes in ACE data, fission reactions actually have
// // an angular distribution listed, but for those that do, it's simply just
// // a uniform distribution in mu
// mu = 2.0 * prn() - 1.0
void sample_fission_neutron(int i_nuclide, const Reaction* rx, double E_in, Bank* site)
{
// Sample cosine of angle -- fission neutrons are always emitted
// isotropically. Sometimes in ACE data, fission reactions actually have
// an angular distribution listed, but for those that do, it's simply just
// a uniform distribution in mu
double mu = 2.0 * prn() - 1.0;
// // Sample azimuthal angle uniformly in [0,2*pi)
// phi = 2.0*PI*prn()
// site % uvw(1) = mu
// site % uvw(2) = std::sqrt(1.0 - mu*mu) * std::cos(phi)
// site % uvw(3) = std::sqrt(1.0 - mu*mu) * std::sin(phi)
// Sample azimuthal angle uniformly in [0,2*pi)
double phi = 2.0*PI*prn();
site->uvw[0] = mu;
site->uvw[1] = std::sqrt(1.0 - mu*mu) * std::cos(phi);
site->uvw[2] = std::sqrt(1.0 - mu*mu) * std::sin(phi);
// // Determine total nu, delayed nu, and delayed neutron fraction
// nu_t = nuc % nu(E_in, EMISSION_TOTAL)
// nu_d = nuc % nu(E_in, EMISSION_DELAYED)
// beta = nu_d / nu_t
// Determine total nu, delayed nu, and delayed neutron fraction
const auto& nuc {data::nuclides[i_nuclide-1]};
double nu_t = nuc->nu(E_in, Nuclide::EmissionMode::total);
double nu_d = nuc->nu(E_in, Nuclide::EmissionMode::delayed);
double beta = nu_d / nu_t;
// if (prn() < beta) {
// // ====================================================================
// // DELAYED NEUTRON SAMPLED
if (prn() < beta) {
// ====================================================================
// DELAYED NEUTRON SAMPLED
// // sampled delayed precursor group
// xi = prn()*nu_d
// prob = 0.0
// do group = 1, nuc % n_precursor
// sampled delayed precursor group
double xi = prn()*nu_d;
double prob = 0.0;
int group;
for (group = 1; group < nuc->n_precursor_; ++group) {
// determine delayed neutron precursor yield for group j
double yield = (*rx->products_[group].yield_)(E_in);
// // determine delayed neutron precursor yield for group j
// yield = rxn % product_yield(1 + group, E_in)
// Check if this group is sampled
prob += yield;
if (xi < prob) break;
}
// // Check if this group is sampled
// prob = prob + yield
// if (xi < prob) exit
// end do
// if the sum of the probabilities is slightly less than one and the
// random number is greater, j will be greater than nuc %
// n_precursor -- check for this condition
group = std::min(group, nuc->n_precursor_);
// // if the sum of the probabilities is slightly less than one and the
// // random number is greater, j will be greater than nuc %
// // n_precursor -- check for this condition
// group = min(group, nuc % n_precursor)
// set the delayed group for the particle born from fission
site->delayed_group = group;
// // set the delayed group for the particle born from fission
// site % delayed_group = group
int n_sample = 0;
while (true) {
// sample from energy/angle distribution -- note that mu has already been
// sampled above and doesn't need to be resampled
rx->products_[group].sample(E_in, site->E, mu);
// n_sample = 0
// do
// // sample from energy/angle distribution -- note that mu has already been
// // sampled above and doesn't need to be resampled
// rxn % product_sample(1 + group, E_in, site % E, mu)
// resample if energy is greater than maximum neutron energy
constexpr int neutron = static_cast<int>(ParticleType::neutron);
if (site->E < data::energy_max[neutron]) break;
// // resample if energy is greater than maximum neutron energy
// if (site % E < energy_max(NEUTRON)) exit
// check for large number of resamples
++n_sample;
if (n_sample == MAX_SAMPLE) {
// particle_write_restart(p)
fatal_error("Resampled energy distribution maximum number of times "
"for nuclide " + nuc->name_);
}
}
// // check for large number of resamples
// n_sample = n_sample + 1
// if (n_sample == MAX_SAMPLE) {
// // particle_write_restart(p)
// fatal_error("Resampled energy distribution maximum number of " &
// // "times for nuclide " // nuc % name)
// }
// end do
} else {
// ====================================================================
// PROMPT NEUTRON SAMPLED
// } else {
// // ====================================================================
// // PROMPT NEUTRON SAMPLED
// set the delayed group for the particle born from fission to 0
site->delayed_group = 0;
// // set the delayed group for the particle born from fission to 0
// site % delayed_group = 0
// sample from prompt neutron energy distribution
int n_sample = 0;
while (true) {
rx->products_[0].sample(E_in, site->E, mu);
// // sample from prompt neutron energy distribution
// n_sample = 0
// do
// rxn % product_sample(1, E_in, site % E, mu)
// resample if energy is greater than maximum neutron energy
constexpr int neutron = static_cast<int>(ParticleType::neutron);
if (site->E < data::energy_max[neutron]) break;
// // resample if energy is greater than maximum neutron energy
// if (site % E < energy_max(NEUTRON)) exit
// // check for large number of resamples
// n_sample = n_sample + 1
// if (n_sample == MAX_SAMPLE) {
// // particle_write_restart(p)
// fatal_error("Resampled energy distribution maximum number of " &
// // "times for nuclide " // nuc % name)
// }
// end do
// }
// }
// check for large number of resamples
++n_sample;
if (n_sample == MAX_SAMPLE) {
// particle_write_restart(p)
fatal_error("Resampled energy distribution maximum number of times "
"for nuclide " + nuc->name_);
}
}
}
}
// void inelastic_scatter(int i_nuclide, const Reaction& rx, Particle* p)
// {