From 1f90ec1560fb2dcf82ce82f493549c8934b5f265 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Mon, 19 Nov 2018 07:47:27 -0600 Subject: [PATCH] Convert sample_fission and sample_fission_neutron to C++ --- include/openmc/nuclide.h | 14 +++ include/openmc/physics.h | 6 +- src/nuclide.cpp | 95 +++++++++++++++- src/physics.F90 | 180 ----------------------------- src/physics.cpp | 239 +++++++++++++++++++-------------------- 5 files changed, 222 insertions(+), 312 deletions(-) diff --git a/include/openmc/nuclide.h b/include/openmc/nuclide.h index b524d4bf7f..5291259ad9 100644 --- a/include/openmc/nuclide.h +++ b/include/openmc/nuclide.h @@ -11,7 +11,9 @@ #include #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 kTs_; //! temperatures in eV (k*T) + bool fissionable_ {false}; //! Whether nuclide is fissionable + bool has_partial_fission_ {false}; //! has partial fission reactions? + std::vector fission_rx_; //! Fission reactions + int n_precursor_ {0}; //! Number of delayed neutron precursors + std::unique_ptr total_nu_; //! Total neutron yield + std::vector> reactions_; //! Reactions private: diff --git a/include/openmc/physics.h b/include/openmc/physics.h index 6ea043c36e..8d1b15978f 100644 --- a/include/openmc/physics.h +++ b/include/openmc/physics.h @@ -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); diff --git a/src/nuclide.cpp b/src/nuclide.cpp index 18a6a29200..4c8cfd313d 100644 --- a/src/nuclide.cpp +++ b/src/nuclide.cpp @@ -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(nu_dset); + } else if (func_type == "Polynomial") { + total_nu_ = std::make_unique(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); } } } diff --git a/src/physics.F90 b/src/physics.F90 index 7f956be907..aaf27ab2e2 100644 --- a/src/physics.F90 +++ b/src/physics.F90 @@ -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. diff --git a/src/physics.cpp b/src/physics.cpp index 44b141b0dc..d19cfdb70a 100644 --- a/src/physics.cpp +++ b/src/physics.cpp @@ -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(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(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) // {