From 88597145922d2974a092776f09957b02649c8c24 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Sat, 17 Nov 2018 15:43:10 -0600 Subject: [PATCH] Move sample_neutron_reaction and create_fission_sites to C++ --- include/openmc/nuclide.h | 73 +++++++----- include/openmc/physics.h | 17 ++- src/nuclide.cpp | 32 ++++- src/nuclide_header.F90 | 7 +- src/physics.F90 | 249 +++++---------------------------------- src/physics.cpp | 190 ++++++++++++++++++++--------- src/simulation.cpp | 2 + 7 files changed, 260 insertions(+), 310 deletions(-) diff --git a/include/openmc/nuclide.h b/include/openmc/nuclide.h index bc3c2a1cc7..b524d4bf7f 100644 --- a/include/openmc/nuclide.h +++ b/include/openmc/nuclide.h @@ -30,24 +30,13 @@ public: int A_; //! Mass number int metastable_; //! Metastable state double awr_; //! Atomic weight ratio + bool fissionable_ {false}; //! Whether nuclide is fissionable std::vector> reactions_; //! Reactions + +private: + void create_derived(); }; -//============================================================================== -// Global variables -//============================================================================== - -namespace data { - -// Minimum/maximum transport energy for each particle type. Order corresponds to -// that of the ParticleType enum -extern std::array energy_min; -extern std::array energy_max; - -extern std::vector> nuclides; - -} // namespace data - //=============================================================================== //! Cached microscopic cross sections for a particular nuclide at the current //! energy @@ -91,19 +80,49 @@ struct NuclideMicroXS { // particle is traveling through //=============================================================================== - struct MaterialMacroXS { - double total; //!< macroscopic total xs - double absorption; //!< macroscopic absorption xs - double fission; //!< macroscopic fission xs - double nu_fission; //!< macroscopic production xs - double photon_prod; //!< macroscopic photon production xs +struct MaterialMacroXS { + double total; //!< macroscopic total xs + double absorption; //!< macroscopic absorption xs + double fission; //!< macroscopic fission xs + double nu_fission; //!< macroscopic production xs + double photon_prod; //!< macroscopic photon production xs - // Photon cross sections - double coherent; //!< macroscopic coherent xs - double incoherent; //!< macroscopic incoherent xs - double photoelectric; //!< macroscopic photoelectric xs - double pair_production; //!< macroscopic pair production xs - }; + // Photon cross sections + double coherent; //!< macroscopic coherent xs + double incoherent; //!< macroscopic incoherent xs + double photoelectric; //!< macroscopic photoelectric xs + double pair_production; //!< macroscopic pair production xs +}; + +//============================================================================== +// Global variables +//============================================================================== + +namespace data { + +// Minimum/maximum transport energy for each particle type. Order corresponds to +// that of the ParticleType enum +extern std::array energy_min; +extern std::array energy_max; + +extern std::vector> nuclides; + +} // namespace data + +namespace simulation { + +// Cross section caches +extern NuclideMicroXS* micro_xs; +extern "C" MaterialMacroXS material_xs; +#pragma omp threadprivate(micro_xs, material_xs) + +} // namespace simulation + +//============================================================================== +// Fortran compatibility +//============================================================================== + +void set_micro_xs(); } // namespace openmc diff --git a/include/openmc/physics.h b/include/openmc/physics.h index 4705065172..dda8b4d7e8 100644 --- a/include/openmc/physics.h +++ b/include/openmc/physics.h @@ -36,17 +36,22 @@ extern "C" void sample_electron_reaction(Particle* p); //! MeV) are created and travel in opposite directions. extern "C" void sample_positron_reaction(Particle* p); -// void sample_nuclide(Particle* p, int mt, int i_nuclide, int i_nuc_mat); +extern "C" 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. +extern "C" void create_fission_sites(Particle* p, int i_nuclide, int i_rx, + Bank* bank_array, int64_t* bank_size, int64_t bank_capacity); // void sample_element(Particle* p); -// int sample_fission(int i_nuclide, double E); +extern "C" int sample_fission(int i_nuclide, double E); // void sample_photon_product(int i_nuclide, double E, int* i_rx, int* i_product); -// void absorption(Particle* p, int i_nuclide); +extern "C" void absorption(Particle* p, int i_nuclide); -// void scatter(Particle*, int i_nuclide, int i_nuc_mat); +extern "C" void scatter(Particle*, int i_nuclide, int i_nuc_mat); // void elastic_scatter(int i_nuclide, const Reaction& rx, double kT, double* E, // Direction* u, double* mu_lab, double* wgt); @@ -59,11 +64,11 @@ extern "C" void sample_positron_reaction(Particle* p); // void sample_cxs_target_velocity(int i_nuclide, Direction* v_target, double E, Direction u, // double kT); -// void sample_fission_neutron(int i_nuclide, const Reaction& rx, double E_in, Bank* site); +extern "C" void sample_fission_neutron(int i_nuclide, int i_rx, double E_in, Bank* site); // void inelastic_scatter(int i_nuclide, const Reaction& rx, Particle* p); -// void sample_secondary_photons(Particle* p, int i_nuclide); +extern "C" void sample_secondary_photons(Particle* p, int i_nuclide); } // namespace openmc diff --git a/src/nuclide.cpp b/src/nuclide.cpp index 629cd86f56..18a6a29200 100644 --- a/src/nuclide.cpp +++ b/src/nuclide.cpp @@ -1,6 +1,7 @@ #include "openmc/nuclide.h" #include "openmc/container_util.h" +#include "openmc/endf.h" #include "openmc/error.h" #include "openmc/hdf5_interface.h" #include "openmc/message_passing.h" @@ -17,14 +18,16 @@ namespace openmc { //============================================================================== namespace data { - std::array energy_min {0.0, 0.0}; std::array energy_max {INFTY, INFTY}; - std::vector> nuclides; - } // namespace data +namespace simulation { +NuclideMicroXS* micro_xs; +MaterialMacroXS material_xs; +} // namespace simulation + //============================================================================== // Nuclide implementation //============================================================================== @@ -152,6 +155,19 @@ Nuclide::Nuclide(hid_t group, const double* temperature, int n) } close_group(rxs_group); + this->create_derived(); +} + +void Nuclide::create_derived() +{ + for (const auto& rx : reactions_) { + // Skip redundant reactions + if (rx->redundant_) continue; + + if (is_fission(rx->mt_)) { + fissionable_ = true; + } + } } //============================================================================== @@ -178,4 +194,14 @@ extern "C" Reaction* nuclide_reaction(Nuclide* nuc, int i_rx) extern "C" void nuclides_clear() { data::nuclides.clear(); } + +extern "C" NuclideMicroXS* micro_xs_ptr(); +void set_micro_xs() +{ +#pragma omp parallel + { + simulation::micro_xs = micro_xs_ptr(); + } +} + } // namespace openmc diff --git a/src/nuclide_header.F90 b/src/nuclide_header.F90 index 3fb19b60f2..7cbe9b85c2 100644 --- a/src/nuclide_header.F90 +++ b/src/nuclide_header.F90 @@ -186,7 +186,7 @@ module nuclide_header ! Cross section caches type(NuclideMicroXS), allocatable, target :: micro_xs(:) ! Cache for each nuclide - type(MaterialMacroXS) :: material_xs ! Cache for current material + type(MaterialMacroXS), bind(C) :: material_xs ! Cache for current material !$omp threadprivate(micro_xs, material_xs) ! Minimum/maximum energies @@ -233,6 +233,11 @@ contains b = library_present_c(type, to_c_string(name)) end function + function micro_xs_ptr() result(ptr) bind(C) + type(C_PTR) :: ptr + ptr = C_LOC(micro_xs(1)) + end function + !=============================================================================== ! ASSIGN_0K_ELASTIC_SCATTERING !=============================================================================== diff --git a/src/physics.F90 b/src/physics.F90 index fd6a3e8ced..045c0d37f5 100644 --- a/src/physics.F90 +++ b/src/physics.F90 @@ -34,91 +34,6 @@ module physics contains -!=============================================================================== -! SAMPLE_NEUTRON_REACTION samples a nuclide based on the macroscopic cross -! sections for each nuclide within a material and then samples a reaction for -! that nuclide and calls the appropriate routine to process the physics. Note -! that there is special logic when suvival biasing is turned on since fission -! and disappearance are treated implicitly. -!=============================================================================== - - subroutine sample_neutron_reaction(p) bind(C) - - type(Particle), intent(inout) :: p - - integer :: i_nuclide ! index in nuclides array - integer :: i_nuc_mat ! index in material's nuclides array - integer :: i_reaction ! index in nuc % reactions array - integer(C_INT) :: err - integer(C_INT64_T) :: n - type(Nuclide), pointer :: nuc - type(C_PTR) :: ptr - type(Bank), pointer :: fission_bank(:) - - call sample_nuclide(p, 'total ', i_nuclide, i_nuc_mat) - - ! Get pointer to table - nuc => nuclides(i_nuclide) - - ! Save which nuclide particle had collision with - p % event_nuclide = i_nuclide - - ! Create fission bank sites. Note that while a fission reaction is sampled, - ! it never actually "happens", i.e. the weight of the particle does not - ! change when sampling fission sites. The following block handles all - ! absorption (including fission) - - if (nuc % fissionable) then - if (run_mode == MODE_EIGENVALUE) then - ! Get fission bank pointer - err = openmc_fission_bank(ptr, n) - call c_f_pointer(ptr, fission_bank, [n]) - - call sample_fission(i_nuclide, p % E, i_reaction) - call create_fission_sites(p, i_nuclide, i_reaction, fission_bank, n_bank) - elseif (run_mode == MODE_FIXEDSOURCE .and. create_fission_neutrons) then - call sample_fission(i_nuclide, p % E, i_reaction) - call create_fission_sites(p, i_nuclide, i_reaction, & - p % secondary_bank, p % n_secondary) - end if - end if - - ! Create secondary photons - if (photon_transport) then - call prn_set_stream(STREAM_PHOTON) - call sample_secondary_photons(p, i_nuclide) - call prn_set_stream(STREAM_TRACKING) - end if - - ! If survival biasing is being used, the following subroutine adjusts the - ! weight of the particle. Otherwise, it checks to see if absorption occurs - - if (micro_xs(i_nuclide) % absorption > ZERO) then - call absorption(p, i_nuclide) - else - p % absorb_wgt = ZERO - end if - if (.not. p % alive) return - - ! Sample a scattering reaction and determine the secondary energy of the - ! exiting neutron - call scatter(p, i_nuclide, i_nuc_mat) - - ! Advance URR seed stream 'N' times after energy changes - if (p % E /= p % last_E) then - call prn_set_stream(STREAM_URR_PTABLE) - call advance_prn_seed(size(nuclides, kind=8)) - call prn_set_stream(STREAM_TRACKING) - end if - - ! Play russian roulette if survival biasing is turned on - if (survival_biasing) then - call russian_roulette(p) - if (.not. p % alive) return - end if - - end subroutine sample_neutron_reaction - !=============================================================================== ! SAMPLE_PHOTON_REACTION samples an element based on the macroscopic cross ! sections for each nuclide within a material and then samples a reaction for @@ -362,12 +277,12 @@ contains ! SAMPLE_NUCLIDE !=============================================================================== - subroutine sample_nuclide(p, base, i_nuclide, i_nuc_mat) + subroutine sample_nuclide(p, base, i_nuclide, i_nuc_mat) bind(C) type(Particle), intent(in) :: p - character(7), intent(in) :: base ! which reaction to sample based on - integer, intent(out) :: i_nuclide - integer, intent(out) :: i_nuc_mat + integer(C_INT), value :: base ! which reaction to sample based on + integer(C_INT), intent(out) :: i_nuclide + integer(C_INT), intent(out) :: i_nuc_mat real(8) :: prob real(8) :: cutoff @@ -380,11 +295,11 @@ contains ! Sample cumulative distribution function select case (base) - case ('total') + case (SCORE_TOTAL) cutoff = prn() * material_xs % total - case ('scatter') + case (SCORE_SCATTER) cutoff = prn() * (material_xs % total - material_xs % absorption) - case ('fission') + case (SCORE_FISSION) cutoff = prn() * material_xs % fission end select @@ -405,12 +320,12 @@ contains ! Determine microscopic cross section select case (base) - case ('total') + case (SCORE_TOTAL) sigma = atom_density * micro_xs(i_nuclide) % total - case ('scatter') + case (SCORE_SCATTER) sigma = atom_density * (micro_xs(i_nuclide) % total - & micro_xs(i_nuclide) % absorption) - case ('fission') + case (SCORE_FISSION) sigma = atom_density * micro_xs(i_nuclide) % fission end select @@ -467,10 +382,10 @@ contains ! SAMPLE_FISSION !=============================================================================== - subroutine sample_fission(i_nuclide, E, i_reaction) - integer, intent(in) :: i_nuclide ! index in nuclides array - real(8), intent(in) :: E ! incident neutron energy - integer, intent(out) :: i_reaction ! index in nuc % reactions array + 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 @@ -529,7 +444,7 @@ contains if (prob > cutoff) exit FISSION_REACTION_LOOP end do FISSION_REACTION_LOOP - end subroutine sample_fission + end function sample_fission !=============================================================================== ! SAMPLE_PHOTON_PRODUCT @@ -594,9 +509,9 @@ contains ! ABSORPTION !=============================================================================== - subroutine absorption(p, i_nuclide) + subroutine absorption(p, i_nuclide) bind(C) type(Particle), intent(inout) :: p - integer, intent(in) :: i_nuclide + integer(C_INT), value :: i_nuclide if (survival_biasing) then ! Determine weight absorbed in survival biasing @@ -634,10 +549,10 @@ contains ! SCATTER !=============================================================================== - subroutine scatter(p, i_nuclide, i_nuc_mat) + subroutine scatter(p, i_nuclide, i_nuc_mat) bind(C) type(Particle), intent(inout) :: p - integer, intent(in) :: i_nuclide - integer, intent(in) :: i_nuc_mat + integer(C_INT), value :: i_nuclide + integer(C_INT), value :: i_nuc_mat integer :: i integer :: j @@ -1137,122 +1052,14 @@ contains end subroutine sample_cxs_target_velocity -!=============================================================================== -! CREATE_FISSION_SITES determines the average total, prompt, and delayed -! neutrons produced from fission and creates appropriate bank sites. -!=============================================================================== - - subroutine create_fission_sites(p, i_nuclide, i_reaction, bank_array, size_bank) - type(Particle), intent(inout) :: p - integer, intent(in) :: i_nuclide - integer, intent(in) :: i_reaction - type(Bank), intent(inout) :: bank_array(:) - integer(8), intent(inout) :: size_bank - - integer :: nu_d(MAX_DELAYED_GROUPS) ! number of delayed neutrons born - integer :: i ! loop index - integer :: nu ! actual number of neutrons produced - real(8) :: nu_t ! total nu - real(8) :: weight ! weight adjustment for ufs method - type(Nuclide), pointer :: nuc - - interface - function ufs_get_weight(p) result(weight) bind(C) - import Particle, C_DOUBLE - type(Particle), intent(in) :: p - real(C_DOUBLE) :: WEIGHT - end function - end interface - - ! Get pointers - nuc => nuclides(i_nuclide) - - ! TODO: Heat generation from fission - - ! If uniform fission source weighting is turned on, we increase of decrease - ! the expected number of fission sites produced - - if (ufs) then - weight = ufs_get_weight(p) - else - weight = ONE - end if - - ! Determine expected number of neutrons produced - nu_t = p % wgt / keff * weight * micro_xs(i_nuclide) % nu_fission / & - micro_xs(i_nuclide) % total - - ! Sample number of neutrons produced - if (prn() > nu_t - int(nu_t)) then - nu = int(nu_t) - else - nu = int(nu_t) + 1 - end if - - ! Check for bank size getting hit. For fixed source calculations, this is a - ! fatal error. For eigenvalue calculations, it just means that k-effective - ! was too high for a single batch. - if (size_bank + nu > size(bank_array)) then - if (run_mode == MODE_FIXEDSOURCE) then - call fatal_error("Secondary particle bank size limit reached. If you & - &are running a subcritical multiplication problem, k-effective & - &may be too close to one.") - else - if (master) call warning("Maximum number of sites in fission bank & - &reached. This can result in irreproducible results using different & - &numbers of processes/threads.") - end if - end if - - ! Bank source neutrons - if (nu == 0 .or. size_bank == size(bank_array)) return - - ! Initialize counter of delayed neutrons encountered for each delayed group - ! to zero. - nu_d(:) = 0 - - p % fission = .true. ! Fission neutrons will be banked - do i = int(size_bank,4) + 1, int(min(size_bank + nu, int(size(bank_array),8)),4) - ! Bank source neutrons by copying particle data - bank_array(i) % xyz = p % coord(1) % xyz - - ! Set particle as neutron - bank_array(i) % particle = NEUTRON - - ! Set weight of fission bank site - bank_array(i) % wgt = ONE/weight - - ! Sample delayed group and angle/energy for fission reaction - call sample_fission_neutron(nuc, nuc % reactions(i_reaction), & - p % E, bank_array(i)) - - ! Set delayed group on particle too - p % delayed_group = bank_array(i) % delayed_group - - ! Increment the number of neutrons born delayed - if (p % delayed_group > 0) then - nu_d(p % delayed_group) = nu_d(p % delayed_group) + 1 - end if - end do - - ! increment number of bank sites - size_bank = min(size_bank + nu, int(size(bank_array),8)) - - ! Store total and delayed weight banked for analog fission tallies - p % n_bank = nu - p % wgt_bank = nu/weight - p % n_delayed_bank(:) = nu_d(:) - - end subroutine create_fission_sites - !=============================================================================== ! SAMPLE_FISSION_NEUTRON !=============================================================================== - subroutine sample_fission_neutron(nuc, rxn, E_in, site) - type(Nuclide), intent(in) :: nuc - type(Reaction), intent(in) :: rxn - real(8), intent(in) :: E_in + 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 @@ -1266,6 +1073,8 @@ contains 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 @@ -1351,6 +1160,8 @@ contains end do end if + end associate + end subroutine sample_fission_neutron !=============================================================================== @@ -1422,9 +1233,9 @@ contains ! SAMPLE_SECONDARY_PHOTONS !=============================================================================== - subroutine sample_secondary_photons(p, i_nuclide) + subroutine sample_secondary_photons(p, i_nuclide) bind(C) type(Particle), intent(inout) :: p - integer, intent(in) :: i_nuclide + integer(C_INT), value :: i_nuclide integer :: i_reaction ! index in nuc % reactions array integer :: i_product ! index in nuc % reactions % products array diff --git a/src/physics.cpp b/src/physics.cpp index b01e40b752..76642858e4 100644 --- a/src/physics.cpp +++ b/src/physics.cpp @@ -1,14 +1,19 @@ #include "openmc/physics.h" +#include "openmc/bank.h" +#include "openmc/constants.h" +#include "openmc/eigenvalue.h" #include "openmc/error.h" +#include "openmc/message_passing.h" #include "openmc/nuclide.h" #include "openmc/photon.h" +#include "openmc/physics_common.h" #include "openmc/random_lcg.h" #include "openmc/reaction.h" #include "openmc/settings.h" #include "openmc/simulation.h" -#include // for max +#include // for max, min #include // for sqrt, exp, log #include @@ -60,70 +65,147 @@ void collision(Particle* p) } } -// void sample_neutron_reaction(Particle* p) -// { -// sample_nuclide(p, 'total ', i_nuclide, i_nuc_mat) +void sample_neutron_reaction(Particle* p) +{ + int i_nuclide; + int i_nuc_mat; + sample_nuclide(p, SCORE_TOTAL, &i_nuclide, &i_nuc_mat); -// // Get pointer to table -// nuc => nuclides(i_nuclide) + // Save which nuclide particle had collision with + p->event_nuclide = i_nuclide; -// // Save which nuclide particle had collision with -// p->event_nuclide = i_nuclide + // Create fission bank sites. Note that while a fission reaction is sampled, + // it never actually "happens", i.e. the weight of the particle does not + // change when sampling fission sites. The following block handles all + // absorption (including fission) -// // Create fission bank sites. Note that while a fission reaction is sampled, -// // it never actually "happens", i.e. the weight of the particle does not -// // change when sampling fission sites. The following block handles all -// // absorption (including fission) + const auto& nuc {data::nuclides[i_nuclide-1]}; -// if (nuc % fissionable) { -// if (run_mode == MODE_EIGENVALUE) { -// // Get fission bank pointer -// err = openmc_fission_bank(ptr, n) -// c_f_pointer(ptr, fission_bank, [n]) + if (nuc->fissionable_) { + int i_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(), + &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, + &p->n_secondary, MAX_SECONDARY); + } + } -// sample_fission(i_nuclide, p->E, i_reaction) -// create_fission_sites(p, i_nuclide, i_reaction, fission_bank, n_bank) -// } else if (run_mode == MODE_FIXEDSOURCE && create_fission_neutrons) { -// sample_fission(i_nuclide, p->E, i_reaction) -// create_fission_sites(p, i_nuclide, i_reaction, & -// p->secondary_bank, p->n_secondary) -// } -// } + // Create secondary photons + if (settings::photon_transport) { + prn_set_stream(STREAM_PHOTON); + sample_secondary_photons(p, i_nuclide); + prn_set_stream(STREAM_TRACKING); + } -// // Create secondary photons -// if (photon_transport) { -// prn_set_stream(STREAM_PHOTON) -// sample_secondary_photons(p, i_nuclide) -// prn_set_stream(STREAM_TRACKING) -// } + // If survival biasing is being used, the following subroutine adjusts the + // weight of the particle. Otherwise, it checks to see if absorption occurs -// // If survival biasing is being used, the following subroutine adjusts the -// // weight of the particle. Otherwise, it checks to see if absorption occurs + if (simulation::micro_xs[i_nuclide-1].absorption > 0.0) { + absorption(p, i_nuclide); + } else { + p->absorb_wgt = 0.0; + } + if (!p->alive) return; -// if (micro_xs(i_nuclide) % absorption > 0.0) { -// absorption(p, i_nuclide) -// } else { -// p->absorb_wgt = 0.0 -// } -// if (!p->alive) return + // Sample a scattering reaction and determine the secondary energy of the + // exiting neutron + scatter(p, i_nuclide, i_nuc_mat); -// // Sample a scattering reaction and determine the secondary energy of the -// // exiting neutron -// scatter(p, i_nuclide, i_nuc_mat) + // Advance URR seed stream 'N' times after energy changes + if (p->E != p->last_E) { + prn_set_stream(STREAM_URR_PTABLE); + advance_prn_seed(data::nuclides.size()); + prn_set_stream(STREAM_TRACKING); + } -// // Advance URR seed stream 'N' times after energy changes -// if (p->E /= p->last_E) { -// prn_set_stream(STREAM_URR_PTABLE) -// advance_prn_seed(size(nuclides, kind=8)) -// prn_set_stream(STREAM_TRACKING) -// } + // Play russian roulette if survival biasing is turned on + if (settings::survival_biasing) { + russian_roulette(p); + if (!p->alive) return; + } +} -// // Play russian roulette if survival biasing is turned on -// if (survival_biasing) { -// russian_roulette(p) -// if (!p->alive) return -// } -// } +void +create_fission_sites(Particle* p, int i_nuclide, int i_rx, Bank* bank_array, + int64_t* size_bank, int64_t bank_capacity) +{ + // TODO: Heat generation from fission + + // If uniform fission source weighting is turned on, we increase or decrease + // the expected number of fission sites produced + double weight = settings::ufs_on ? ufs_get_weight(p) : 1.0; + + // Determine the expected number of neutrons produced + double nu_t = p->wgt / simulation::keff * weight * simulation::micro_xs[ + i_nuclide-1].nu_fission / simulation::micro_xs[i_nuclide-1].total; + + // Sample the number of neutrons produced + int nu = static_cast(nu_t); + if (prn() <= (nu_t - nu)) ++nu; + + // Check for the bank size getting hit. For fixed source calculations, this + // is a fatal error; for eigenvalue calculations, it just means that k-eff + // was too high for a single batch. + if (*size_bank + nu > bank_capacity) { + if (settings::run_mode == RUN_MODE_FIXEDSOURCE) { + throw std::runtime_error{"Secondary particle bank size limit reached." + " If you are running a subcritical multiplication problem," + " k-effective may be too close to one."}; + } else { + if (mpi::master) { + warning("Maximum number of sites in fission bank reached. This can" + " result in irreproducible results using different numbers of" + " processes/threads."); + } + } + } + + // Begin banking the source neutrons + // First, if our bank is full then don't continue + if (nu == 0 || *size_bank == bank_capacity) return; + + // Initialize the counter of delayed neutrons encountered for each delayed + // group. + double nu_d[MAX_DELAYED_GROUPS] = {0.}; + + p->fission = true; + for (size_t i = *size_bank; i < std::min(*size_bank + nu, bank_capacity); ++i) { + // Bank source neutrons by copying the particle data + bank_array[i].xyz[0] = p->coord[0].xyz[0]; + bank_array[i].xyz[1] = p->coord[0].xyz[1]; + bank_array[i].xyz[2] = p->coord[0].xyz[2]; + + // Set that the bank particle is a neutron + bank_array[i].particle = static_cast(ParticleType::neutron); + + // Set the weight of the fission bank site + 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]); + + // Set the delayed group on the particle as well + p->delayed_group = bank_array[i].delayed_group; + + // Increment the number of neutrons born delayed + if (p->delayed_group > 0) { + nu_d[p->delayed_group-1]++; + } + } + + // Increment number of bank sites + *size_bank = std::min(*size_bank + nu, bank_capacity); + + // Store the total weight banked for analog fission tallies + p->n_bank = nu; + p->wgt_bank = nu / weight; + for (size_t d = 0; d < MAX_DELAYED_GROUPS; d++) { + p->n_delayed_bank[d] = nu_d[d]; + } +} // void sample_photon_reaction(Particle* p) // { diff --git a/src/simulation.cpp b/src/simulation.cpp index 0cfb99a4ae..8459c9c71f 100644 --- a/src/simulation.cpp +++ b/src/simulation.cpp @@ -6,6 +6,7 @@ #include "openmc/eigenvalue.h" #include "openmc/error.h" #include "openmc/message_passing.h" +#include "openmc/nuclide.h" #include "openmc/output.h" #include "openmc/particle.h" #include "openmc/random_lcg.h" @@ -100,6 +101,7 @@ int openmc_simulation_init() // Call Fortran initialization simulation_init_f(); + set_micro_xs(); // Reset global variables -- this is done before loading state point (as that // will potentially populate k_generation and entropy)