From e8af7a2d16e8ca219992561fe7ad3770ece3bb2a Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Wed, 21 Nov 2018 11:13:55 -0600 Subject: [PATCH] Convert scattering routines to C++ --- include/openmc/endf.h | 11 + include/openmc/material.h | 6 + include/openmc/nuclide.h | 26 +- include/openmc/physics.h | 30 +- include/openmc/position.h | 6 + src/endf.cpp | 29 ++ src/input_xml.F90 | 4 - src/material_header.F90 | 17 + src/nuclide.cpp | 111 +++++- src/nuclide_header.F90 | 162 --------- src/physics.F90 | 567 +----------------------------- src/physics.cpp | 721 ++++++++++++++++++++------------------ src/position.cpp | 18 + src/settings.cpp | 12 - 14 files changed, 620 insertions(+), 1100 deletions(-) diff --git a/include/openmc/endf.h b/include/openmc/endf.h index c781bd0b84..a7030d0355 100644 --- a/include/openmc/endf.h +++ b/include/openmc/endf.h @@ -22,6 +22,17 @@ Interpolation int2interp(int i); //! \return Whether corresponding reaction is a fission reaction bool is_fission(int MT); +//! Determine if a given MT number is that of a disappearance reaction, i.e., a +//! reaction with no neutron in the exit channel +//! \param[in] MT ENDF MT value +//! \return Whether corresponding reaction is a disappearance reaction +bool is_disappearance(int MT); + +//! Determine if a given MT number is that of an inelastic scattering reaction +//! \param[in] MT ENDF MT value +//! \return Whether corresponding reaction is an inelastic scattering reaction +bool is_inelastic_scatter(int MT); + //============================================================================== //! Abstract one-dimensional function //============================================================================== diff --git a/include/openmc/material.h b/include/openmc/material.h index fcab023a52..ad9b9c9f5a 100644 --- a/include/openmc/material.h +++ b/include/openmc/material.h @@ -43,5 +43,11 @@ public: explicit Material(pugi::xml_node material_node); }; +//============================================================================== +// Fortran compatibility +//============================================================================== + +extern "C" bool material_isotropic(int i_material, int i_nuc_mat); + } // namespace openmc #endif // OPENMC_MATERIAL_H diff --git a/include/openmc/nuclide.h b/include/openmc/nuclide.h index 5291259ad9..348dae2c75 100644 --- a/include/openmc/nuclide.h +++ b/include/openmc/nuclide.h @@ -17,19 +17,35 @@ namespace openmc { +//============================================================================== +// Constants +//============================================================================== + +constexpr double CACHE_INVALID {-1.0}; + //=============================================================================== // Data for a nuclide //=============================================================================== class Nuclide { public: + // Types, aliases using EmissionMode = ReactionProduct::EmissionMode; + struct EnergyGrid { + std::vector grid_index; + std::vector energy; + }; // Constructors Nuclide(hid_t group, const double* temperature, int n); // Methods - double nu(double E, EmissionMode mode, int group=0); + double nu(double E, EmissionMode mode, int group=0) const; + void calculate_elastic_xs(int i_nuclide) const; + + //! Determines the microscopic 0K elastic cross section at a trial relative + //! energy used in resonance scattering + double elastic_xs_0K(double E) const; // Data members std::string name_; //! Name of nuclide, e.g. "U235" @@ -38,6 +54,7 @@ public: int metastable_; //! Metastable state double awr_; //! Atomic weight ratio std::vector kTs_; //! temperatures in eV (k*T) + std::vector grid_; //! Energy grid at each temperature bool fissionable_ {false}; //! Whether nuclide is fissionable bool has_partial_fission_ {false}; //! has partial fission reactions? @@ -45,7 +62,14 @@ public: int n_precursor_ {0}; //! Number of delayed neutron precursors std::unique_ptr total_nu_; //! Total neutron yield + // Resonance scattering information + bool resonant_ {false}; + std::vector energy_0K_; + std::vector elastic_0K_; + std::vector xs_cdf_; + std::vector> reactions_; //! Reactions + std::vector index_inelastic_scatter_; private: void create_derived(); diff --git a/include/openmc/physics.h b/include/openmc/physics.h index b428bf8238..58bc7e9906 100644 --- a/include/openmc/physics.h +++ b/include/openmc/physics.h @@ -2,6 +2,7 @@ #define OPENMC_PHYSICS_H #include "openmc/bank.h" +#include "openmc/nuclide.h" #include "openmc/particle.h" #include "openmc/position.h" #include "openmc/reaction.h" @@ -51,22 +52,33 @@ 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 scatter(Particle*, int i_nuclide, int i_nuc_mat); +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); +//! Treats the elastic scattering of a neutron with a target. +void elastic_scatter(int i_nuclide, const Reaction* rx, double kT, double* E, + double* uvw, double* mu_lab, double* wgt); -// void sab_scatter(int i_nuclide, int i_sab, double* E, Direction* u, double* mu); +extern "C" void sab_scatter(int i_nuclide, int i_sab, double* E, + double* uvw, double* mu); -// void sample_target_velocity(int i_nuclide, Direction* v_target, double E, Direction u, -// Direction v_neut, double* wgt, double xs_eff, double kT); +//! samples the target velocity. The constant cross section free gas model is +//! the default method. Methods for correctly accounting for the energy +//! dependence of cross sections in treating resonance elastic scattering such +//! as the DBRC, WCM, and a new, accelerated scheme are also implemented here. +Direction sample_target_velocity(const Nuclide* nuc, double E, Direction u, + Direction v_neut, double xs_eff, double kT, double* wgt); -// void sample_cxs_target_velocity(int i_nuclide, Direction* v_target, double E, Direction u, -// double kT); +//! samples a target velocity based on the free gas scattering formulation, used +//! by most Monte Carlo codes, in which cross section is assumed to be constant +//! in energy. Excellent documentation for this method can be found in +//! FRA-TM-123. +Direction sample_cxs_target_velocity(double awr, double E, Direction u, double kT); 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); +//! handles all reactions with a single secondary neutron (other than fission), +//! i.e. level scattering, (n,np), (n,na), etc. +void inelastic_scatter(const Nuclide* nuc, const Reaction* rx, Particle* p); void sample_secondary_photons(Particle* p, int i_nuclide); diff --git a/include/openmc/position.h b/include/openmc/position.h index 62c7f804fc..0deb298dd2 100644 --- a/include/openmc/position.h +++ b/include/openmc/position.h @@ -25,6 +25,8 @@ struct Position { Position& operator-=(double); Position& operator*=(Position); Position& operator*=(double); + Position& operator/=(Position); + Position& operator/=(double); const double& operator[](int i) const { switch (i) { @@ -76,6 +78,10 @@ inline Position operator*(Position a, Position b) { return a *= b; } inline Position operator*(Position a, double b) { return a *= b; } inline Position operator*(double a, Position b) { return b *= a; } +inline Position operator/(Position a, Position b) { return a /= b; } +inline Position operator/(Position a, double b) { return a /= b; } +inline Position operator/(double a, Position b) { return b /= a; } + inline bool operator==(Position a, Position b) {return a.x == b.x && a.y == b.y && a.z == b.z;} diff --git a/src/endf.cpp b/src/endf.cpp index 8ffb75345d..284f14c08f 100644 --- a/src/endf.cpp +++ b/src/endf.cpp @@ -45,6 +45,35 @@ bool is_fission(int mt) return mt == 18 || mt == 19 || mt == 20 || mt == 21 || mt == 38; } +bool is_disappearance(int mt) +{ + if (mt >= N_DISAPPEAR && mt <= N_DA) { + return true; + } else if (mt >= N_P0 && mt <= N_AC) { + return true; + } else if (mt == N_TA || mt == N_DT || mt == N_P3HE || mt == N_D3HE + || mt == N_3HEA || mt == N_3P) { + return true; + } else { + return false; + } +} + +bool is_inelastic_scatter(int mt) +{ + if (mt < 100) { + if (is_fission(mt)) { + return false; + } else { + return mt >= MISC && mt != 27; + } + } else if (mt <= 200) { + return !is_disappearance(mt); + } else { + return false; + } +} + //============================================================================== // Polynomial implementation //============================================================================== diff --git a/src/input_xml.F90 b/src/input_xml.F90 index 1da1bf9df6..dab02f52bc 100644 --- a/src/input_xml.F90 +++ b/src/input_xml.F90 @@ -2113,9 +2113,6 @@ contains call close_group(group_id) call file_close(file_id) - ! Assign resonant scattering data - if (res_scat_on) call nuclides(i_nuclide) % assign_0K_elastic_scattering() - ! Determine if minimum/maximum energy for this nuclide is greater/less ! than the previous if (size(nuclides(i_nuclide) % grid) >= 1) then @@ -2297,7 +2294,6 @@ contains logical :: file_exists ! Does multipole library exist? character(7) :: readable ! Is multipole library readable? character(MAX_FILE_LEN) :: filename ! Path to multipole xs library - character(kind=C_CHAR), pointer :: string(:) integer(HID_T) :: file_id integer(HID_T) :: group_id diff --git a/src/material_header.F90 b/src/material_header.F90 index 3bcf96db8b..e041f52f4a 100644 --- a/src/material_header.F90 +++ b/src/material_header.F90 @@ -1041,4 +1041,21 @@ contains end subroutine bremsstrahlung_init +!=============================================================================== +! Fortran compatibility +!=============================================================================== + + function material_isotropic(i_material, i_nuc_mat) result(iso) bind(C) + integer(C_INT), value :: i_material + integer(C_INT), value :: i_nuc_mat + logical(C_BOOL) :: iso + + iso = .false. + associate (mat => materials(i_material)) + if (mat % has_isotropic_nuclides) then + iso = mat % p0(i_nuc_mat) + end if + end associate + end function + end module material_header diff --git a/src/nuclide.cpp b/src/nuclide.cpp index 4c8cfd313d..77b3177cd8 100644 --- a/src/nuclide.cpp +++ b/src/nuclide.cpp @@ -5,6 +5,7 @@ #include "openmc/error.h" #include "openmc/hdf5_interface.h" #include "openmc/message_passing.h" +#include "openmc/search.h" #include "openmc/settings.h" #include "openmc/string_utils.h" @@ -143,22 +144,48 @@ 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 + hid_t energy_group = open_group(group, "energy"); for (const auto& T : temps_to_read) { std::string dset {std::to_string(T) + "K"}; + + // Determine exact kT values double kT; read_dataset(kT_group, dset.c_str(), kT); kTs_.push_back(kT); + + // Read energy grid + grid_.emplace_back(); + read_dataset(energy_group, dset.c_str(), grid_.back().energy); } close_group(kT_group); + // Check for 0K energy grid + if (object_exists(energy_group, "0K")) { + read_dataset(energy_group, "0K", energy_0K_); + } + close_group(energy_group); + // Read reactions hid_t rxs_group = open_group(group, "reactions"); for (auto name : group_names(rxs_group)) { if (starts_with(name, "reaction_")) { hid_t rx_group = open_group(rxs_group, name.c_str()); reactions_.push_back(std::make_unique(rx_group, temps_to_read)); + + // Check for 0K elastic scattering + if (reactions_.back()->mt_ == ELASTIC) { + if (object_exists(rx_group, "0K")) { + hid_t temp_group = open_group(rx_group, "0K"); + read_dataset(temp_group, "xs", elastic_0K_); + close_group(temp_group); + } + } close_group(rx_group); + + // Determine reaction indices for inelastic scattering reactions + if (is_inelastic_scatter(reactions_.back()->mt_)) { + index_inelastic_scatter_.push_back(reactions_.size() - 1); + } } } close_group(rxs_group); @@ -211,9 +238,52 @@ void Nuclide::create_derived() } } } + + if (settings::res_scat_on) { + // Determine if this nuclide should be treated as a resonant scatterer + if (!settings::res_scat_nuclides.empty()) { + // If resonant nuclides were specified, check the list explicitly + for (const auto& name : settings::res_scat_nuclides) { + if (name_ == name) { + resonant_ = true; + + // Make sure nuclide has 0K data + if (energy_0K_.empty()) { + fatal_error("Cannot treat " + name_ + " as a resonant scatterer " + "because 0 K elastic scattering data is not present."); + } + break; + } + } + } else { + // Otherwise, assume that any that have 0 K elastic scattering data are + // resonant + resonant_ = !energy_0K_.empty(); + } + + if (resonant_) { + // Build CDF for 0K elastic scattering + double xs_cdf_sum = 0.0; + xs_cdf_.resize(energy_0K_.size()); + xs_cdf_[0] = 0.0; + + const auto& E = energy_0K_; + auto& xs = elastic_0K_; + for (int i = 0; i < E.size() - 1; ++i) { + // Negative cross sections result in a CDF that is not monotonically + // increasing. Set all negative xs values to zero. + if (xs[i] < 0.0) xs[i] = 0.0; + + // build xs cdf + xs_cdf_sum += (std::sqrt(E[i])*xs[i] + std::sqrt(E[i+1])*xs[i+1]) + / 2.0 * (E[i+1] - E[i]); + xs_cdf_[i] = xs_cdf_sum; + } + } + } } -double Nuclide::nu(double E, EmissionMode mode, int group) +double Nuclide::nu(double E, EmissionMode mode, int group) const { if (!fissionable_) return 0.0; @@ -253,6 +323,43 @@ double Nuclide::nu(double E, EmissionMode mode, int group) } } +void Nuclide::calculate_elastic_xs(int i_nuclide) const +{ + // Get temperature index, grid index, and interpolation factor + auto& micro = simulation::micro_xs[i_nuclide-1]; + int i_temp = micro.index_temp - 1; + int i_grid = micro.index_grid - 1; + double f = micro.interp_factor; + + if (i_temp >= 0) { + const auto& xs = reactions_[0]->xs_[i_temp].value; + micro.elastic = (1.0 - f)*xs[i_grid] + f*xs[i_grid + 1]; + } +} + +double Nuclide::elastic_xs_0K(double E) const +{ + // Determine index on nuclide energy grid + int i_grid; + if (E < energy_0K_.front()) { + i_grid = 0; + } else if (E > energy_0K_.back()) { + i_grid = energy_0K_.size() - 2; + } else { + i_grid = lower_bound_index(energy_0K_.begin(), energy_0K_.end(), E); + } + + // check for rare case where two energy points are the same + if (energy_0K_[i_grid] == energy_0K_[i_grid+1]) ++i_grid; + + // calculate interpolation factor + double f = (E - energy_0K_[i_grid]) / + (energy_0K_[i_grid + 1] - energy_0K_[i_grid]); + + // Calculate microscopic nuclide elastic cross section + return (1.0 - f)*elastic_0K_[i_grid] + f*elastic_0K_[i_grid + 1]; +} + //============================================================================== // Fortran compatibility functions //============================================================================== diff --git a/src/nuclide_header.F90 b/src/nuclide_header.F90 index 3427be8cf7..5db1355db2 100644 --- a/src/nuclide_header.F90 +++ b/src/nuclide_header.F90 @@ -71,15 +71,7 @@ module nuclide_header ! Microscopic cross sections type(SumXS), allocatable :: xs(:) - ! Resonance scattering info - logical :: resonant = .false. ! resonant scatterer? - real(8), allocatable :: energy_0K(:) ! energy grid for 0K xs - real(8), allocatable :: elastic_0K(:) ! Microscopic elastic cross section - real(8), allocatable :: xs_cdf(:) ! CDF of v_rel times cross section - ! Fission information - logical :: has_partial_fission = .false. ! nuclide has partial fission reactions? - integer :: n_fission = 0 ! # of fission reactions integer :: n_precursor = 0 ! # of delayed neutron precursors integer, allocatable :: index_fission(:) ! indices in reactions class(Function1D), allocatable :: total_nu @@ -95,7 +87,6 @@ module nuclide_header ! Reactions type(Reaction), allocatable :: reactions(:) - integer, allocatable :: index_inelastic_scatter(:) ! Array that maps MT values to index in reactions; used at tally-time. Note ! that ENDF-102 does not assign any MT values above 891. @@ -108,7 +99,6 @@ module nuclide_header type(C_PTR) :: ptr contains - procedure :: assign_0K_elastic_scattering procedure :: clear => nuclide_clear procedure :: from_hdf5 => nuclide_from_hdf5 procedure :: init_grid => nuclide_init_grid @@ -238,79 +228,6 @@ contains ptr = C_LOC(micro_xs(1)) end function -!=============================================================================== -! ASSIGN_0K_ELASTIC_SCATTERING -!=============================================================================== - - subroutine assign_0K_elastic_scattering(this) - class(Nuclide), intent(inout) :: this - - integer :: i - real(8) :: xs_cdf_sum - - interface - function res_scat_nuclides_empty() result(empty) bind(C) - import C_BOOL - logical(C_BOOL) :: empty - end function - - function res_scat_nuclides_size() result(n) bind(C) - import C_INT - integer(C_INT) :: n - end function - - function res_scat_nuclides_cmp(i, name) result(b) bind(C) - import C_INT, C_CHAR, C_BOOL - integer(C_INT), value :: i - character(kind=C_CHAR), intent(in) :: name(*) - logical(C_BOOL) :: b - end function - end interface - - this % resonant = .false. - if (.not. res_scat_nuclides_empty()) then - ! If resonant nuclides were specified, check the list explicitly - do i = 1, res_scat_nuclides_size() - if (res_scat_nuclides_cmp(i, to_c_string(this % name))) then - this % resonant = .true. - - ! Make sure nuclide has 0K data - if (.not. allocated(this % energy_0K)) then - call fatal_error("Cannot treat " // trim(this % name) // " as a & - &resonant scatterer because 0 K elastic scattering data is & - ¬ present.") - end if - - exit - end if - end do - else - ! Otherwise, assume that any that have 0 K elastic scattering data are - ! resonant - this % resonant = allocated(this % energy_0K) - end if - - if (this % resonant) then - ! Build CDF for 0K elastic scattering - xs_cdf_sum = ZERO - allocate(this % xs_cdf(0:size(this % energy_0K))) - this % xs_cdf(0) = ZERO - - associate (E => this % energy_0K, xs => this % elastic_0K) - do i = 1, size(E) - 1 - ! Negative cross sections result in a CDF that is not monotonically - ! increasing. Set all negative xs values to zero. - if (xs(i) < ZERO) xs(i) = ZERO - - ! build xs cdf - xs_cdf_sum = xs_cdf_sum + (sqrt(E(i))*xs(i) + sqrt(E(i+1))*xs(i+1))& - / TWO * (E(i+1) - E(i)) - this % xs_cdf(i) = xs_cdf_sum - end do - end associate - end if - end subroutine assign_0K_elastic_scattering - !=============================================================================== ! NUCLIDE_CLEAR resets and deallocates data in Nuclide !=============================================================================== @@ -341,7 +258,6 @@ contains integer(HID_T) :: kT_group integer(HID_T) :: rxs_group integer(HID_T) :: rx_group - integer(HID_T) :: xs, temp_group integer(HID_T) :: total_nu integer(HID_T) :: fer_group ! fission_energy_release group integer(HID_T) :: fer_dset @@ -355,7 +271,6 @@ contains real(8) :: temp_actual type(VectorInt) :: MTs type(VectorInt) :: temps_to_read - type(VectorInt) :: index_inelastic_scatter interface function nuclide_from_hdf5_c(group, temperature, n) result(ptr) bind(C) @@ -494,15 +409,6 @@ contains call read_dataset(this % grid(i) % energy, energy_dset) call close_dataset(energy_dset) end do - - ! Check for 0K energy grid - if (object_exists(energy_group, '0K')) then - energy_dset = open_dataset(energy_group, '0K') - call get_shape(energy_dset, dims) - allocate(this % energy_0K(int(dims(1), 4))) - call read_dataset(this % energy_0K, energy_dset) - call close_dataset(energy_dset) - end if call close_group(energy_group) ! Get MT values based on group names @@ -523,34 +429,10 @@ contains ! Set pointer for each reaction call this % reactions(i) % init(this % ptr, i) - ! Check for 0K elastic scattering - if (this % reactions(i) % MT == 2) then - if (object_exists(rx_group, '0K')) then - temp_group = open_group(rx_group, '0K') - xs = open_dataset(temp_group, 'xs') - call get_shape(xs, dims) - allocate(this % elastic_0K(int(dims(1), 4))) - call read_dataset(this % elastic_0K, xs) - call close_dataset(xs) - call close_group(temp_group) - end if - end if - - ! Add the reaction index to the scattering array if this is an inelastic - ! scatter reaction - if (is_inelastic_scatter(MTs % data(i))) then - call index_inelastic_scatter % push_back(i) - end if - call close_group(rx_group) end do call close_group(rxs_group) - ! Recast to a regular array to save space - allocate(this % index_inelastic_scatter(index_inelastic_scatter % size())) - this % index_inelastic_scatter = & - index_inelastic_scatter % data(1: index_inelastic_scatter % size()) - ! Read unresolved resonance probability tables if present if (object_exists(group_id, 'urr')) then this % urr_present = .true. @@ -726,7 +608,6 @@ contains allocate(this % index_fission(1)) elseif (rx % MT == N_F) then allocate(this % index_fission(PARTIAL_FISSION_MAX)) - this % has_partial_fission = .true. end if end if @@ -746,7 +627,6 @@ contains if (t == 1) then i_fission = i_fission + 1 this % index_fission(i_fission) = i - this % n_fission = this % n_fission + 1 end if end if ! fission end do ! temperature @@ -1361,45 +1241,6 @@ contains end if end subroutine multipole_deriv_eval -!=============================================================================== -! 0K_ELASTIC_XS determines the microscopic 0K elastic cross section -! for a given nuclide at the trial relative energy used in resonance scattering -!=============================================================================== - - pure function elastic_xs_0K(E, nuc) result(xs_out) - real(8), intent(in) :: E ! trial energy - type(Nuclide), intent(in) :: nuc ! target nuclide at temperature - real(8) :: xs_out ! 0K xs at trial energy - - integer :: i_grid ! index on nuclide energy grid - integer :: n_grid - real(8) :: f ! interp factor on nuclide energy grid - - ! Determine index on nuclide energy grid - n_grid = size(nuc % energy_0K) - if (E < nuc % energy_0K(1)) then - i_grid = 1 - elseif (E > nuc % energy_0K(n_grid)) then - i_grid = n_grid - 1 - else - i_grid = binary_search(nuc % energy_0K, n_grid, E) - end if - - ! check for rare case where two energy points are the same - if (nuc % energy_0K(i_grid) == nuc % energy_0K(i_grid+1)) then - i_grid = i_grid + 1 - end if - - ! calculate interpolation factor - f = (E - nuc % energy_0K(i_grid)) & - & / (nuc % energy_0K(i_grid + 1) - nuc % energy_0K(i_grid)) - - ! Calculate microscopic nuclide elastic cross section - xs_out = (ONE - f) * nuc % elastic_0K(i_grid) & - & + f * nuc % elastic_0K(i_grid + 1) - - end function elastic_xs_0K - !=============================================================================== ! CALCULATE_URR_XS determines cross sections in the unresolved resonance range ! from probability tables @@ -1678,9 +1519,6 @@ contains call nuclide_dict % set(to_lower(name_), n) n_nuclides = n - ! Assign resonant scattering data - if (res_scat_on) call nuclides(n) % assign_0K_elastic_scattering() - ! Initialize nuclide grid call nuclides(n) % init_grid(energy_min(NEUTRON), & energy_max(NEUTRON), n_log_bins) diff --git a/src/physics.F90 b/src/physics.F90 index 8f9acbf903..c01b17142e 100644 --- a/src/physics.F90 +++ b/src/physics.F90 @@ -1,9 +1,6 @@ module physics - use algorithm, only: binary_search - use bank_header use constants - use endf, only: reaction_name use error, only: fatal_error, warning, write_message use material_header, only: Material, materials use math @@ -15,12 +12,10 @@ module physics atomic_relaxation, pair_production, & thick_target_bremsstrahlung use physics_common - use random_lcg, only: prn, advance_prn_seed, prn_set_stream - use reaction_header, only: Reaction + use random_lcg, only: prn use sab_header, only: sab_tables use settings use simulation_header - use string, only: to_str use tally_header implicit none @@ -254,229 +249,17 @@ contains end function sample_element -!=============================================================================== -! SCATTER -!=============================================================================== - - subroutine scatter(p, i_nuclide, i_nuc_mat) bind(C) - type(Particle), intent(inout) :: p - integer(C_INT), value :: i_nuclide - integer(C_INT), value :: i_nuc_mat - - integer :: i - integer :: j - integer :: i_temp - integer :: i_grid - integer :: threshold - real(8) :: f - real(8) :: prob - real(8) :: cutoff - real(8) :: uvw_new(3) ! outgoing uvw for iso-in-lab scattering - real(8) :: uvw_old(3) ! incoming uvw for iso-in-lab scattering - real(8) :: phi ! azimuthal angle for iso-in-lab scattering - real(8) :: kT ! temperature in eV - logical :: sampled ! whether or not a reaction type has been sampled - type(Nuclide), pointer :: nuc - - ! copy incoming direction - uvw_old(:) = p % coord(1) % uvw - - ! Get pointer to nuclide and grid index/interpolation factor - nuc => nuclides(i_nuclide) - i_temp = micro_xs(i_nuclide) % index_temp - i_grid = micro_xs(i_nuclide) % index_grid - f = micro_xs(i_nuclide) % interp_factor - - ! For tallying purposes, this routine might be called directly. In that - ! case, we need to sample a reaction via the cutoff variable - cutoff = prn() * (micro_xs(i_nuclide) % total - & - micro_xs(i_nuclide) % absorption) - sampled = .false. - - ! Calculate elastic cross section if it wasn't precalculated - if (micro_xs(i_nuclide) % elastic == CACHE_INVALID) then - call nuc % calculate_elastic_xs(micro_xs(i_nuclide)) - end if - - prob = micro_xs(i_nuclide) % elastic - micro_xs(i_nuclide) % thermal - if (prob > cutoff) then - ! ======================================================================= - ! NON-S(A,B) ELASTIC SCATTERING - - ! Determine temperature - if (nuc % mp_present) then - kT = p % sqrtkT**2 - else - kT = nuc % kTs(micro_xs(i_nuclide) % index_temp) - end if - - ! Perform collision physics for elastic scattering - call elastic_scatter(i_nuclide, nuc % reactions(1), kT, p % E, & - p % coord(1) % uvw, p % mu, p % wgt) - - p % event_MT = ELASTIC - sampled = .true. - end if - - prob = micro_xs(i_nuclide) % elastic - if (prob > cutoff .and. .not. sampled) then - ! ======================================================================= - ! S(A,B) SCATTERING - - call sab_scatter(i_nuclide, micro_xs(i_nuclide) % index_sab, p % E, & - p % coord(1) % uvw, p % mu) - - p % event_MT = ELASTIC - sampled = .true. - end if - - if (.not. sampled) then - ! ======================================================================= - ! INELASTIC SCATTERING - - j = 0 - do while (prob < cutoff) - j = j + 1 - i = nuc % index_inelastic_scatter(j) - - ! Check to make sure inelastic scattering reaction sampled - if (i > size(nuc % reactions)) then - call particle_write_restart(p) - call fatal_error("Did not sample any reaction for nuclide " & - &// trim(nuc % name)) - end if - - associate (rx => nuc % reactions(i)) - ! 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 - end do - - ! Perform collision physics for inelastic scattering - call inelastic_scatter(nuc, nuc%reactions(i), p) - p % event_MT = nuc % reactions(i) % MT - - end if - - ! Set event component - p % event = EVENT_SCATTER - - ! Sample new outgoing angle for isotropic-in-lab scattering - associate (mat => materials(p % material)) - if (mat % has_isotropic_nuclides) then - if (materials(p % material) % p0(i_nuc_mat)) then - ! Sample isotropic-in-lab outgoing direction - uvw_new(1) = TWO * prn() - ONE - phi = TWO * PI * prn() - uvw_new(2) = cos(phi) * sqrt(ONE - uvw_new(1)*uvw_new(1)) - uvw_new(3) = sin(phi) * sqrt(ONE - uvw_new(1)*uvw_new(1)) - p % mu = dot_product(uvw_old, uvw_new) - - ! Change direction of particle - p % coord(1) % uvw = uvw_new - end if - end if - end associate - - end subroutine scatter - -!=============================================================================== -! ELASTIC_SCATTER treats the elastic scattering of a neutron with a -! target. -!=============================================================================== - - subroutine elastic_scatter(i_nuclide, rxn, kT, E, uvw, mu_lab, wgt) - integer, intent(in) :: i_nuclide - type(Reaction), intent(in) :: rxn - real(8), intent(in) :: kT ! temperature in eV - real(8), intent(inout) :: E - real(8), intent(inout) :: uvw(3) - real(8), intent(out) :: mu_lab - real(8), intent(inout) :: wgt - - real(8) :: awr ! atomic weight ratio of target - real(8) :: mu_cm ! cosine of polar angle in center-of-mass - real(8) :: vel ! magnitude of velocity - real(8) :: v_n(3) ! velocity of neutron - real(8) :: v_cm(3) ! velocity of center-of-mass - real(8) :: v_t(3) ! velocity of target nucleus - real(8) :: uvw_cm(3) ! directional cosines in center-of-mass - type(Nuclide), pointer :: nuc - - ! get pointer to nuclide - nuc => nuclides(i_nuclide) - - vel = sqrt(E) - awr = nuc % awr - - ! Neutron velocity in LAB - v_n = vel * uvw - - ! Sample velocity of target nucleus - if (.not. micro_xs(i_nuclide) % use_ptable) then - call sample_target_velocity(nuc, v_t, E, uvw, v_n, wgt, & - micro_xs(i_nuclide) % elastic, kT) - else - v_t = ZERO - end if - - ! Velocity of center-of-mass - v_cm = (v_n + awr*v_t)/(awr + ONE) - - ! Transform to CM frame - v_n = v_n - v_cm - - ! Find speed of neutron in CM - vel = sqrt(dot_product(v_n, v_n)) - - ! Sample scattering angle - mu_cm = rxn % sample_elastic_mu(E) - - ! Determine direction cosines in CM - uvw_cm = v_n/vel - - ! Rotate neutron velocity vector to new angle -- note that the speed of the - ! neutron in CM does not change in elastic scattering. However, the speed - ! will change when we convert back to LAB - v_n = vel * rotate_angle(uvw_cm, mu_cm) - - ! Transform back to LAB frame - v_n = v_n + v_cm - - E = dot_product(v_n, v_n) - vel = sqrt(E) - - ! compute cosine of scattering angle in LAB frame by taking dot product of - ! neutron's pre- and post-collision angle - mu_lab = dot_product(uvw, v_n) / vel - - ! Set energy and direction of particle in LAB frame - uvw = v_n / vel - - ! Because of floating-point roundoff, it may be possible for mu_lab to be - ! outside of the range [-1,1). In these cases, we just set mu_lab to exactly - ! -1 or 1 - - if (abs(mu_lab) > ONE) mu_lab = sign(ONE,mu_lab) - - end subroutine elastic_scatter - !=============================================================================== ! SAB_SCATTER performs thermal scattering of a particle with a bound scatterer ! according to a specified S(a,b) table. !=============================================================================== - subroutine sab_scatter(i_nuclide, i_sab, E, uvw, mu) - integer, intent(in) :: i_nuclide ! index in micro_xs - integer, intent(in) :: i_sab ! index in sab_tables - real(8), intent(inout) :: E ! incoming/outgoing energy - real(8), intent(inout) :: uvw(3) ! directional cosines - real(8), intent(out) :: mu ! scattering cosine + subroutine sab_scatter(i_nuclide, i_sab, E, uvw, mu) bind(C) + integer(C_INT), value :: i_nuclide ! index in micro_xs + integer(C_INT), value :: i_sab ! index in sab_tables + real(C_DOUBLE), intent(inout) :: E ! incoming/outgoing energy + real(C_DOUBLE), intent(inout) :: uvw(3) ! directional cosines + real(C_DOUBLE), intent(out) :: mu ! scattering cosine real(C_DOUBLE) :: E_out type(C_PTR) :: ptr @@ -490,340 +273,4 @@ contains uvw = rotate_angle(uvw, mu) end subroutine sab_scatter -!=============================================================================== -! SAMPLE_TARGET_VELOCITY samples the target velocity. The constant cross section -! free gas model is the default method. Methods for correctly accounting -! for the energy dependence of cross sections in treating resonance elastic -! scattering such as the DBRC, WCM, and a new, accelerated scheme are also -! implemented here. -!=============================================================================== - - subroutine sample_target_velocity(nuc, v_target, E, uvw, v_neut, wgt, xs_eff, kT) - type(Nuclide), intent(in) :: nuc ! target nuclide at temperature T - real(8), intent(out) :: v_target(3) ! target velocity - real(8), intent(in) :: E ! particle energy - real(8), intent(in) :: uvw(3) ! direction cosines - real(8), intent(in) :: v_neut(3) ! neutron velocity - real(8), intent(inout) :: wgt ! particle weight - real(8), intent(in) :: xs_eff ! effective elastic xs at temperature T - real(8), intent(in) :: kT ! equilibrium temperature of target in eV - - real(8) :: awr ! target/neutron mass ratio - real(8) :: E_rel ! trial relative energy - real(8) :: xs_0K ! 0K xs at E_rel - real(8) :: wcf ! weight correction factor - real(8) :: E_red ! reduced energy (same as used by Cullen in SIGMA1) - real(8) :: E_low ! lowest practical relative energy - real(8) :: E_up ! highest practical relative energy - real(8) :: E_t ! trial target energy - real(8) :: xs_max ! max 0K xs over practical relative energies - real(8) :: xs_low ! 0K xs at lowest practical relative energy - real(8) :: xs_up ! 0K xs at highest practical relative energy - real(8) :: m ! slope for interpolation - real(8) :: R ! rejection criterion for DBRC / target speed - real(8) :: cdf_low ! xs cdf at lowest practical relative energy - real(8) :: cdf_up ! xs cdf at highest practical relative energy - real(8) :: cdf_rel ! trial xs cdf value - real(8) :: mu ! cosine between neutron and target velocities - - integer :: i_E_low ! 0K index to lowest practical relative energy - integer :: i_E_up ! 0K index to highest practical relative energy - integer :: i_E_rel ! index to trial relative energy - integer :: n_grid ! number of energies on 0K grid - - integer :: sampling_method ! method of target velocity sampling - - awr = nuc % awr - - ! check if nuclide is a resonant scatterer - if (nuc % resonant) then - - ! sampling method to use - sampling_method = res_scat_method - - ! upper resonance scattering energy bound (target is at rest above this E) - if (E > res_scat_energy_max) then - v_target = ZERO - return - - ! lower resonance scattering energy bound (should be no resonances below) - else if (E < res_scat_energy_min) then - sampling_method = RES_SCAT_CXS - end if - - ! otherwise, use free gas model - else - if (E >= FREE_GAS_THRESHOLD * kT .and. awr > ONE) then - v_target = ZERO - return - else - sampling_method = RES_SCAT_CXS - end if - end if - - ! use appropriate target velocity sampling method - select case (sampling_method) - case (RES_SCAT_CXS) - - ! sample target velocity with the constant cross section (cxs) approx. - call sample_cxs_target_velocity(nuc, v_target, E, uvw, kT) - - case (RES_SCAT_WCM) - - ! sample target velocity with the constant cross section (cxs) approx. - call sample_cxs_target_velocity(nuc, v_target, E, uvw, kT) - - ! adjust weight as prescribed by the weight correction method (wcm) - E_rel = dot_product((v_neut - v_target), (v_neut - v_target)) - xs_0K = elastic_xs_0K(E_rel, nuc) - wcf = xs_0K / xs_eff - wgt = wcf * wgt - - case (RES_SCAT_DBRC, RES_SCAT_ARES) - E_red = sqrt(awr * E / kT) - E_low = max(ZERO, E_red - FOUR)**2 * kT / awr - E_up = (E_red + FOUR)**2 * kT / awr - - ! find lower and upper energy bound indices - ! lower index - n_grid = size(nuc % energy_0K) - if (E_low < nuc % energy_0K(1)) then - i_E_low = 1 - elseif (E_low > nuc % energy_0K(n_grid)) then - i_E_low = n_grid - 1 - else - i_E_low = binary_search(nuc % energy_0K, n_grid, E_low) - end if - - ! upper index - if (E_up < nuc % energy_0K(1)) then - i_E_up = 1 - elseif (E_up > nuc % energy_0K(n_grid)) then - i_E_up = n_grid - 1 - else - i_E_up = binary_search(nuc % energy_0K, n_grid, E_up) - end if - - if (i_E_up == i_E_low) then - ! Handle degenerate case -- if the upper/lower bounds occur for the same - ! index, then using cxs is probably a good approximation - call sample_cxs_target_velocity(nuc, v_target, E, uvw, kT) - - else - if (sampling_method == RES_SCAT_DBRC) then - ! interpolate xs since we're not exactly at the energy indices - xs_low = nuc % elastic_0K(i_E_low) - m = (nuc % elastic_0K(i_E_low + 1) - xs_low) & - / (nuc % energy_0K(i_E_low + 1) - nuc % energy_0K(i_E_low)) - xs_low = xs_low + m * (E_low - nuc % energy_0K(i_E_low)) - xs_up = nuc % elastic_0K(i_E_up) - m = (nuc % elastic_0K(i_E_up + 1) - xs_up) & - / (nuc % energy_0K(i_E_up + 1) - nuc % energy_0K(i_E_up)) - xs_up = xs_up + m * (E_up - nuc % energy_0K(i_E_up)) - - ! get max 0K xs value over range of practical relative energies - xs_max = max(xs_low, & - maxval(nuc % elastic_0K(i_E_low + 1 : i_E_up)), xs_up) - - DBRC_REJECT_LOOP: do - TARGET_ENERGY_LOOP: do - ! sample target velocity with the constant cross section (cxs) approx. - call sample_cxs_target_velocity(nuc, v_target, E, uvw, kT) - E_rel = dot_product((v_neut - v_target), (v_neut - v_target)) - if (E_rel < E_up) exit TARGET_ENERGY_LOOP - end do TARGET_ENERGY_LOOP - - ! perform Doppler broadening rejection correction (dbrc) - xs_0K = elastic_xs_0K(E_rel, nuc) - R = xs_0K / xs_max - if (prn() < R) exit DBRC_REJECT_LOOP - end do DBRC_REJECT_LOOP - - elseif (sampling_method == RES_SCAT_ARES) then - ! interpolate xs CDF since we're not exactly at the energy indices - ! cdf value at lower bound attainable energy - m = (nuc % xs_cdf(i_E_low) - nuc % xs_cdf(i_E_low - 1)) & - / (nuc % energy_0K(i_E_low + 1) - nuc % energy_0K(i_E_low)) - cdf_low = nuc % xs_cdf(i_E_low - 1) & - + m * (E_low - nuc % energy_0K(i_E_low)) - if (E_low <= nuc % energy_0K(1)) cdf_low = ZERO - - ! cdf value at upper bound attainable energy - m = (nuc % xs_cdf(i_E_up) - nuc % xs_cdf(i_E_up - 1)) & - / (nuc % energy_0K(i_E_up + 1) - nuc % energy_0K(i_E_up)) - cdf_up = nuc % xs_cdf(i_E_up - 1) & - + m * (E_up - nuc % energy_0K(i_E_up)) - - ARES_REJECT_LOOP: do - - ! directly sample Maxwellian - E_t = -kT * log(prn()) - - ! sample a relative energy using the xs cdf - cdf_rel = cdf_low + prn() * (cdf_up - cdf_low) - i_E_rel = binary_search(nuc % xs_cdf(i_E_low-1:i_E_up), & - i_E_up - i_E_low + 2, cdf_rel) - E_rel = nuc % energy_0K(i_E_low + i_E_rel - 1) - m = (nuc % xs_cdf(i_E_low + i_E_rel - 1) & - - nuc % xs_cdf(i_E_low + i_E_rel - 2)) & - / (nuc % energy_0K(i_E_low + i_E_rel) & - - nuc % energy_0K(i_E_low + i_E_rel - 1)) - E_rel = E_rel + (cdf_rel - nuc % xs_cdf(i_E_low + i_E_rel - 2)) / m - - ! perform rejection sampling on cosine between - ! neutron and target velocities - mu = (E_t + awr * (E - E_rel)) / (TWO * sqrt(awr * E * E_t)) - - if (abs(mu) < ONE) then - ! set and accept target velocity - E_t = E_t / awr - v_target = sqrt(E_t) * rotate_angle(uvw, mu) - exit ARES_REJECT_LOOP - end if - end do ARES_REJECT_LOOP - end if - end if - end select - - end subroutine sample_target_velocity - -!=============================================================================== -! SAMPLE_CXS_TARGET_VELOCITY samples a target velocity based on the free gas -! scattering formulation, used by most Monte Carlo codes, in which cross section -! is assumed to be constant in energy. Excellent documentation for this method -! can be found in FRA-TM-123. -!=============================================================================== - - subroutine sample_cxs_target_velocity(nuc, v_target, E, uvw, kT) - type(Nuclide), intent(in) :: nuc ! target nuclide at temperature - real(8), intent(out) :: v_target(3) - real(8), intent(in) :: E - real(8), intent(in) :: uvw(3) - real(8), intent(in) :: kT ! equilibrium temperature of target in eV - - real(8) :: awr ! target/neutron mass ratio - real(8) :: alpha ! probability of sampling f2 over f1 - real(8) :: mu ! cosine of angle between neutron and target vel - real(8) :: r1, r2 ! pseudo-random numbers - real(8) :: c ! cosine used in maxwell sampling - real(8) :: accept_prob ! probability of accepting combination of vt and mu - real(8) :: beta_vn ! beta * speed of neutron - real(8) :: beta_vt ! beta * speed of target - real(8) :: beta_vt_sq ! (beta * speed of target)^2 - real(8) :: vt ! speed of target - - awr = nuc % awr - - beta_vn = sqrt(awr * E / kT) - alpha = ONE/(ONE + sqrt(pi)*beta_vn/TWO) - - do - ! Sample two random numbers - r1 = prn() - r2 = prn() - - if (prn() < alpha) then - ! With probability alpha, we sample the distribution p(y) = - ! y*e^(-y). This can be done with sampling scheme C45 frmo the Monte - ! Carlo sampler - - beta_vt_sq = -log(r1*r2) - - else - ! With probability 1-alpha, we sample the distribution p(y) = y^2 * - ! e^(-y^2). This can be done with sampling scheme C61 from the Monte - ! Carlo sampler - - c = cos(PI/TWO * prn()) - beta_vt_sq = -log(r1) - log(r2)*c*c - end if - - ! Determine beta * vt - beta_vt = sqrt(beta_vt_sq) - - ! Sample cosine of angle between neutron and target velocity - mu = TWO*prn() - ONE - - ! Determine rejection probability - accept_prob = sqrt(beta_vn*beta_vn + beta_vt_sq - 2*beta_vn*beta_vt*mu) & - /(beta_vn + beta_vt) - - ! Perform rejection sampling on vt and mu - if (prn() < accept_prob) exit - end do - - ! Determine speed of target nucleus - vt = sqrt(beta_vt_sq*kT/awr) - - ! Determine velocity vector of target nucleus based on neutron's velocity - ! and the sampled angle between them - v_target = vt * rotate_angle(uvw, mu) - - end subroutine sample_cxs_target_velocity - -!=============================================================================== -! INELASTIC_SCATTER handles all reactions with a single secondary neutron (other -! than fission), i.e. level scattering, (n,np), (n,na), etc. -!=============================================================================== - - subroutine inelastic_scatter(nuc, rxn, p) - type(Nuclide), intent(in) :: nuc - type(Reaction), intent(in) :: rxn - type(Particle), intent(inout) :: p - - integer :: i ! loop index - real(8) :: E ! energy in lab (incoming/outgoing) - real(8) :: mu ! cosine of scattering angle in lab - real(8) :: A ! atomic weight ratio of nuclide - real(8) :: E_in ! incoming energy - real(8) :: E_cm ! outgoing energy in center-of-mass - real(8) :: yield ! neutron yield - - ! copy energy of neutron - E_in = p % E - - ! sample outgoing energy and scattering cosine - 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 - if (rxn % scatter_in_cm) then - E_cm = E - - ! determine outgoing energy in lab - A = nuc%awr - E = E_cm + (E_in + TWO * mu * (A+ONE) * sqrt(E_in * E_cm)) & - / ((A+ONE)*(A+ONE)) - - ! determine outgoing angle in lab - mu = mu * sqrt(E_cm/E) + ONE/(A+ONE) * sqrt(E_in/E) - end if - - ! Because of floating-point roundoff, it may be possible for mu to be - ! outside of the range [-1,1). In these cases, we just set mu to exactly -1 - ! or 1 - if (abs(mu) > ONE) mu = sign(ONE,mu) - - ! Set outgoing energy and scattering angle - p % E = E - p % mu = mu - - ! change direction of particle - p % coord(1) % uvw = rotate_angle(p % coord(1) % uvw, mu) - - ! evaluate yield - 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 - call particle_create_secondary(p, p % coord(1) % uvw, p % E, & - NEUTRON, run_CE=.true._C_BOOL) - end do - else - ! Otherwise, change weight of particle based on yield - p % wgt = yield * p % wgt - end if - - end subroutine inelastic_scatter - end module physics diff --git a/src/physics.cpp b/src/physics.cpp index 482d92599c..98e215399d 100644 --- a/src/physics.cpp +++ b/src/physics.cpp @@ -4,6 +4,7 @@ #include "openmc/constants.h" #include "openmc/eigenvalue.h" #include "openmc/error.h" +#include "openmc/material.h" #include "openmc/math_functions.h" #include "openmc/message_passing.h" #include "openmc/nuclide.h" @@ -11,12 +12,15 @@ #include "openmc/physics_common.h" #include "openmc/random_lcg.h" #include "openmc/reaction.h" +#include "openmc/secondary_uncorrelated.h" +#include "openmc/search.h" #include "openmc/settings.h" #include "openmc/simulation.h" +#include "openmc/thermal.h" #include "openmc/tallies/tally.h" -#include // for max, min -#include // for sqrt, exp, log +#include // for max, min, max_element +#include // for sqrt, exp, log, abs, copysign #include namespace openmc { @@ -509,24 +513,24 @@ Reaction* sample_fission(int i_nuclide, double E) // } // Get grid index and interpolatoin factor and sample fission cdf - int i_temp = simulation::micro_xs[i_nuclide-1].index_temp; + int i_temp = simulation::micro_xs[i_nuclide-1].index_temp - 1; 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 - for (auto& rx : nuc->reactions_) { + for (auto& rx : nuc->fission_rx_) { // if energy is below threshold for this reaction, skip it - int threshold = rx->xs_[i_temp-1].threshold; + int threshold = rx->xs_[i_temp].threshold; if (i_grid < threshold) continue; // add to cumulative probability - prob += (1.0 - f) * rx->xs_[i_temp-1].value[i_grid - threshold] + prob += (1.0 - f) * rx->xs_[i_temp].value[i_grid - threshold] + f*rx->xs_[i_temp].value[i_grid - threshold + 1]; // Create fission bank sites if fission occurs - if (prob > cutoff) break; + if (prob > cutoff) return rx; } } @@ -596,173 +600,183 @@ void absorption(Particle* p, int i_nuclide) } } -// void scatter(Particle*, int i_nuclide, int i_nuc_mat) -// { -// // copy incoming direction -// uvw_old(:) = p->coord(1) % uvw +void scatter(Particle* p, int i_nuclide, int i_nuc_mat) +{ + // copy incoming direction + Direction u_old {p->coord[0].uvw}; -// // Get pointer to nuclide and grid index/interpolation factor -// nuc => nuclides(i_nuclide) -// i_temp = simulation::micro_xs[i_nuclide-1].index_temp -// i_grid = simulation::micro_xs[i_nuclide-1].index_grid -// f = simulation::micro_xs[i_nuclide-1].interp_factor + // Get pointer to nuclide and grid index/interpolation factor + const auto& nuc {data::nuclides[i_nuclide-1]}; + const auto& micro {simulation::micro_xs[i_nuclide-1]}; + int i_temp = micro.index_temp - 1; + int i_grid = micro.index_grid - 1; + double f = micro.interp_factor; -// // For tallying purposes, this routine might be called directly. In that -// // case, we need to sample a reaction via the cutoff variable -// cutoff = prn() * (micro_xs[i_nuclide-1].total - & -// simulation::micro_xs[i_nuclide-1].absorption) -// sampled = false + // For tallying purposes, this routine might be called directly. In that + // case, we need to sample a reaction via the cutoff variable + double cutoff = prn() * (micro.total - micro.absorption); + bool sampled = false; -// // Calculate elastic cross section if it wasn't precalculated -// if (micro_xs[i_nuclide-1].elastic == CACHE_INVALID) { -// nuc % calculate_elastic_xs(micro_xs(i_nuclide)) -// } + // Calculate elastic cross section if it wasn't precalculated + if (micro.elastic == CACHE_INVALID) { + nuc->calculate_elastic_xs(i_nuclide); + } -// prob = simulation::micro_xs[i_nuclide-1].elastic - simulation::micro_xs[i_nuclide-1].thermal -// if (prob > cutoff) { -// // ======================================================================= -// // NON-S(A,B) ELASTIC SCATTERING + double prob = micro.elastic - micro.thermal; + if (prob > cutoff) { + // ======================================================================= + // NON-S(A,B) ELASTIC SCATTERING -// // Determine temperature -// if (nuc % mp_present) { -// kT = p->sqrtkT**2 -// } else { -// kT = nuc % kTs(micro_xs[i_nuclide-1].index_temp) -// } + // Determine temperature + double kT; + // if (nuc % mp_present) { + // kT = p->sqrtkT**2 + // } else { + kT = nuc->kTs_[i_temp]; + // } -// // Perform collision physics for elastic scattering -// elastic_scatter(i_nuclide, nuc % reactions(1), kT, p->E, & -// p->coord(1) % uvw, p->mu, p->wgt) + // Perform collision physics for elastic scattering + elastic_scatter(i_nuclide, nuc->reactions_[0].get(), kT, + &p->E, p->coord[0].uvw, &p->mu, &p->wgt); -// p->event_MT = ELASTIC -// sampled = true -// } + p->event_MT = ELASTIC; + sampled = true; + } -// prob = simulation::micro_xs[i_nuclide-1].elastic -// if (prob > cutoff && !sampled) { -// // ======================================================================= -// // S(A,B) SCATTERING + prob = micro.elastic; + if (prob > cutoff && !sampled) { + // ======================================================================= + // S(A,B) SCATTERING -// sab_scatter(i_nuclide, simulation::micro_xs[i_nuclide-1].index_sab, p->E, & -// p->coord(1) % uvw, p->mu) + sab_scatter(i_nuclide, micro.index_sab, &p->E, p->coord[0].uvw, &p->mu); -// p->event_MT = ELASTIC -// sampled = true -// } + p->event_MT = ELASTIC; + sampled = true; + } -// if (!sampled) { -// // ======================================================================= -// // INELASTIC SCATTERING + if (!sampled) { + // ======================================================================= + // INELASTIC SCATTERING -// j = 0 -// do while (prob < cutoff) -// j = j + 1 -// i = nuc % index_inelastic_scatter(j) + int j = 0; + int i; + while (prob < cutoff) { + i = nuc->index_inelastic_scatter_[j]; + ++j; -// // Check to make sure inelastic scattering reaction sampled -// if (i > size(nuc % reactions)) { -// particle_write_restart(p) -// fatal_error("Did not sample any reaction for nuclide " & -// &// trim(nuc % name)) -// } + // Check to make sure inelastic scattering reaction sampled + if (i >= nuc->reactions_.size()) { + p->write_restart(); + fatal_error("Did not sample any reaction for nuclide " + nuc->name_); + } -// associate (rx => nuc % reactions(i)) -// // if energy is below threshold for this reaction, skip it -// threshold = rx % xs_threshold(i_temp) -// if (i_grid < threshold) cycle + // if energy is below threshold for this reaction, skip it + const auto& xs {nuc->reactions_[i]->xs_[i_temp]}; + int threshold = xs.threshold - 1; + if (i_grid < threshold) continue; -// // 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 -// end do + // add to cumulative probability + prob += (1.0 - f)*xs.value[i_grid - threshold] + + f*xs.value[i_grid - threshold + 1]; + } -// // Perform collision physics for inelastic scattering -// inelastic_scatter(nuc, nuc%reactions(i), p) -// p->event_MT = nuc % reactions(i) % MT + // Perform collision physics for inelastic scattering + const auto& rx {nuc->reactions_[i]}; + inelastic_scatter(nuc.get(), rx.get(), p); + p->event_MT = rx->mt_; + } -// } + // Set event component + p->event = EVENT_SCATTER; -// // Set event component -// p->event = EVENT_SCATTER + // Sample new outgoing angle for isotropic-in-lab scattering + if (material_isotropic(p->material, i_nuc_mat)) { + // Sample isotropic-in-lab outgoing direction + double mu = 2.0*prn() - 1.0; + double phi = 2.0*PI*prn(); + Direction u_new; + u_new.x = mu; + u_new.y = std::sqrt(1.0 - mu*mu)*std::cos(phi); + u_new.z = std::sqrt(1.0 - mu*mu)*std::sin(phi); -// // Sample new outgoing angle for isotropic-in-lab scattering -// associate (mat => materials(p->material)) -// if (mat % has_isotropic_nuclides) { -// if (materials(p->material) % p0(i_nuc_mat)) { -// // Sample isotropic-in-lab outgoing direction -// uvw_new(1) = 2.0 * prn() - 1.0 -// phi = 2.0 * PI * prn() -// uvw_new(2) = std::cos(phi) * std::sqrt(1.0 - uvw_new(1)*uvw_new(1)) -// uvw_new(3) = std::sin(phi) * std::sqrt(1.0 - uvw_new(1)*uvw_new(1)) -// p->mu = dot_product(uvw_old, uvw_new) + p->mu = u_old.dot(u_new); -// // Change direction of particle -// p->coord(1) % uvw = uvw_new -// } -// } -// end associate -// } + // Change direction of particle + p->coord[0].uvw[0] = u_new.x; + p->coord[0].uvw[1] = u_new.y; + p->coord[0].uvw[2] = u_new.z; + } +} -// void elastic_scatter(int i_nuclide, const Reaction& rx, double kT, double* E, -// Direction& u, double* mu_lab, double* wgt) -// { -// // get pointer to nuclide -// nuc => nuclides(i_nuclide) +void elastic_scatter(int i_nuclide, const Reaction* rx, double kT, double* E, + double* uvw, double* mu_lab, double* wgt) +{ + // get pointer to nuclide + const auto& nuc {data::nuclides[i_nuclide-1]}; -// vel = std::sqrt(E) -// awr = nuc % awr + double vel = std::sqrt(*E); + double awr = nuc->awr_; -// // Neutron velocity in LAB -// v_n = vel * uvw + // Neutron velocity in LAB + Direction u {uvw}; + Direction v_n = vel*u; -// // Sample velocity of target nucleus -// if (!micro_xs[i_nuclide-1].use_ptable) { -// sample_target_velocity(nuc, v_t, E, uvw, v_n, wgt, & -// simulation::micro_xs[i_nuclide-1].elastic, kT) -// } else { -// v_t = 0.0 -// } + // Sample velocity of target nucleus + Direction v_t {}; + if (!simulation::micro_xs[i_nuclide-1].use_ptable) { + v_t = sample_target_velocity(nuc.get(), *E, u, v_n, + simulation::micro_xs[i_nuclide-1].elastic, kT, wgt); + } -// // Velocity of center-of-mass -// v_cm = (v_n + awr*v_t)/(awr + 1.0) + // Velocity of center-of-mass + Direction v_cm = (v_n + awr*v_t)/(awr + 1.0); -// // Transform to CM frame -// v_n = v_n - v_cm + // Transform to CM frame + v_n -= v_cm; -// // Find speed of neutron in CM -// vel = std::sqrt(dot_product(v_n, v_n)) + // Find speed of neutron in CM + vel = v_n.norm(); -// // Sample scattering angle -// mu_cm = rxn % sample_elastic_mu(E) + // Sample scattering angle, checking if it is an ncorrelated angle-energy + // distribution + double mu_cm; + auto& d = rx->products_[0].distribution_[0]; + auto d_ = dynamic_cast(d.get()); + if (d_) { + mu_cm = d_->angle().sample(*E); + } else { + mu_cm = 2.0*prn() - 1.0; + } -// // Determine direction cosines in CM -// uvw_cm = v_n/vel + // Determine direction cosines in CM + Direction u_cm = v_n/vel; -// // Rotate neutron velocity vector to new angle -- note that the speed of the -// // neutron in CM does not change in elastic scattering. However, the speed -// // will change when we convert back to LAB -// v_n = vel * rotate_angle(uvw_cm, mu_cm) + // Rotate neutron velocity vector to new angle -- note that the speed of the + // neutron in CM does not change in elastic scattering. However, the speed + // will change when we convert back to LAB + v_n = vel * rotate_angle(u_cm, mu_cm, nullptr); -// // Transform back to LAB frame -// v_n = v_n + v_cm + // Transform back to LAB frame + v_n += v_cm; -// E = dot_product(v_n, v_n) -// vel = std::sqrt(E) + *E = v_n.dot(v_n); + vel = std::sqrt(*E); -// // compute cosine of scattering angle in LAB frame by taking dot product of -// // neutron's pre- and post-collision angle -// mu_lab = dot_product(uvw, v_n) / vel + // compute cosine of scattering angle in LAB frame by taking dot product of + // neutron's pre- and post-collision angle + *mu_lab = u.dot(v_n) / vel; -// // Set energy and direction of particle in LAB frame -// uvw = v_n / vel + // Set energy and direction of particle in LAB frame + u = v_n / vel; + uvw[0] = u.x; + uvw[1] = u.y; + uvw[2] = u.z; -// // Because of floating-point roundoff, it may be possible for mu_lab to be -// // outside of the range [-1,1). In these cases, we just set mu_lab to exactly -// // -1 or 1 - -// if (abs(mu_lab) > 1.0) mu_lab = sign(1.0,mu_lab) -// } + // Because of floating-point roundoff, it may be possible for mu_lab to be + // outside of the range [-1,1). In these cases, we just set mu_lab to exactly + // -1 or 1 + if (std::abs(*mu_lab) > 1.0) *mu_lab = std::copysign(1.0, *mu_lab); +} // void sab_scatter(int i_nuclide, int i_sab, double* E, Direction* u, double* mu) // { @@ -775,212 +789,217 @@ void absorption(Particle* p, int i_nuclide) // uvw = rotate_angle(uvw, mu) // } -// void sample_target_velocity(int i_nuclide, Direction* v_target, double E, Direction u, -// Direction v_neut, double* wgt, double xs_eff, double kT) -// { -// awr = nuc % awr +Direction sample_target_velocity(const Nuclide* nuc, double E, Direction u, + Direction v_neut, double xs_eff, double kT, double* wgt) +{ + // check if nuclide is a resonant scatterer + int sampling_method; + if (nuc->resonant_) { -// // check if nuclide is a resonant scatterer -// if (nuc % resonant) { + // sampling method to use + sampling_method = settings::res_scat_method; -// // sampling method to use -// sampling_method = res_scat_method + // upper resonance scattering energy bound (target is at rest above this E) + if (E > settings::res_scat_energy_max) { + return {}; -// // upper resonance scattering energy bound (target is at rest above this E) -// if (E > res_scat_energy_max) { -// v_target = 0.0 -// return + // lower resonance scattering energy bound (should be no resonances below) + } else if (E < settings::res_scat_energy_min) { + sampling_method = RES_SCAT_CXS; + } -// // lower resonance scattering energy bound (should be no resonances below) -// } else if (E < res_scat_energy_min) { -// sampling_method = RES_SCAT_CXS -// } + // otherwise, use free gas model + } else { + if (E >= FREE_GAS_THRESHOLD * kT && nuc->awr_ > 1.0) { + return {}; + } else { + sampling_method = RES_SCAT_CXS; + } + } -// // otherwise, use free gas model -// } else { -// if (E >= FREE_GAS_THRESHOLD * kT && awr > 1.0) { -// v_target = 0.0 -// return -// } else { -// sampling_method = RES_SCAT_CXS -// } -// } + // use appropriate target velocity sampling method + switch (sampling_method) { + case RES_SCAT_CXS: -// // use appropriate target velocity sampling method -// select case (sampling_method) -// case (RES_SCAT_CXS) + // sample target velocity with the constant cross section (cxs) approx. + return sample_cxs_target_velocity(nuc->awr_, E, u, kT); -// // sample target velocity with the constant cross section (cxs) approx. -// sample_cxs_target_velocity(nuc, v_target, E, uvw, kT) + case RES_SCAT_WCM: { + // sample target velocity with the constant cross section (cxs) approx. + Direction v_target = sample_cxs_target_velocity(nuc->awr_, E, u, kT); -// case (RES_SCAT_WCM) + // adjust weight as prescribed by the weight correction method (wcm) + Direction v_rel = v_neut - v_target; + double E_rel = v_rel.dot(v_rel); + double xs_0K = nuc->elastic_xs_0K(E_rel); + *wgt *= xs_0K / xs_eff; + return v_target; + } -// // sample target velocity with the constant cross section (cxs) approx. -// sample_cxs_target_velocity(nuc, v_target, E, uvw, kT) + case RES_SCAT_DBRC: + case RES_SCAT_ARES: { + double E_red = std::sqrt(nuc->awr_ * E / kT); + double E_low = std::pow(std::max(0.0, E_red - 4.0), 2) * kT / nuc->awr_; + double E_up = (E_red + 4.0)*(E_red + 4.0) * kT / nuc->awr_; -// // adjust weight as prescribed by the weight correction method (wcm) -// E_rel = dot_product((v_neut - v_target), (v_neut - v_target)) -// xs_0K = elastic_xs_0K(E_rel, nuc) -// wcf = xs_0K / xs_eff -// wgt = wcf * wgt + // find lower and upper energy bound indices + // lower index + int i_E_low; + if (E_low < nuc->energy_0K_.front()) { + i_E_low = 0; + } else if (E_low > nuc->energy_0K_.back()) { + i_E_low = nuc->energy_0K_.size() - 2; + } else { + i_E_low = lower_bound_index(nuc->energy_0K_.begin(), + nuc->energy_0K_.end(), E_low); + } -// case (RES_SCAT_DBRC, RES_SCAT_ARES) -// E_red = std::sqrt(awr * E / kT) -// E_low = std::max(0.0, E_red - FOUR)**2 * kT / awr -// E_up = (E_red + FOUR)**2 * kT / awr + // upper index + int i_E_up; + if (E_up < nuc->energy_0K_.front()) { + i_E_up = 0; + } else if (E_up > nuc->energy_0K_.back()) { + i_E_up = nuc->energy_0K_.size() - 2; + } else { + i_E_up = lower_bound_index(nuc->energy_0K_.begin(), + nuc->energy_0K_.end(), E_up); + } -// // find lower and upper energy bound indices -// // lower index -// n_grid = size(nuc % energy_0K) -// if (E_low < nuc % energy_0K(1)) { -// i_E_low = 1 -// } else if (E_low > nuc % energy_0K(n_grid)) { -// i_E_low = n_grid - 1 -// } else { -// i_E_low = binary_search(nuc % energy_0K, n_grid, E_low) -// } + if (i_E_up == i_E_low) { + // Handle degenerate case -- if the upper/lower bounds occur for the same + // index, then using cxs is probably a good approximation + return sample_cxs_target_velocity(nuc->awr_, E, u, kT); + } -// // upper index -// if (E_up < nuc % energy_0K(1)) { -// i_E_up = 1 -// } else if (E_up > nuc % energy_0K(n_grid)) { -// i_E_up = n_grid - 1 -// } else { -// i_E_up = binary_search(nuc % energy_0K, n_grid, E_up) -// } + if (sampling_method == RES_SCAT_DBRC) { + // interpolate xs since we're not exactly at the energy indices + double xs_low = nuc->elastic_0K_[i_E_low]; + double m = (nuc->elastic_0K_[i_E_low + 1] - xs_low) + / (nuc->energy_0K_[i_E_low + 1] - nuc->energy_0K_[i_E_low]); + xs_low += m * (E_low - nuc->energy_0K_[i_E_low]); + double xs_up = nuc->elastic_0K_[i_E_up]; + m = (nuc->elastic_0K_[i_E_up + 1] - xs_up) + / (nuc->energy_0K_[i_E_up + 1] - nuc->energy_0K_[i_E_up]); + xs_up += m * (E_up - nuc->energy_0K_[i_E_up]); -// if (i_E_up == i_E_low) { -// // Handle degenerate case -- if the upper/lower bounds occur for the same -// // index, then using cxs is probably a good approximation -// sample_cxs_target_velocity(nuc, v_target, E, uvw, kT) + // get max 0K xs value over range of practical relative energies + double xs_max = *std::max_element(&nuc->elastic_0K_[i_E_low + 1], + &nuc->elastic_0K_[i_E_up + 1]); + xs_max = std::max({xs_low, xs_max, xs_up}); -// } else { -// if (sampling_method == RES_SCAT_DBRC) { -// // interpolate xs since we're not exactly at the energy indices -// xs_low = nuc % elastic_0K(i_E_low) -// m = (nuc % elastic_0K(i_E_low + 1) - xs_low) & -// / (nuc % energy_0K(i_E_low + 1) - nuc % energy_0K(i_E_low)) -// xs_low = xs_low + m * (E_low - nuc % energy_0K(i_E_low)) -// xs_up = nuc % elastic_0K(i_E_up) -// m = (nuc % elastic_0K(i_E_up + 1) - xs_up) & -// / (nuc % energy_0K(i_E_up + 1) - nuc % energy_0K(i_E_up)) -// xs_up = xs_up + m * (E_up - nuc % energy_0K(i_E_up)) + while (true) { + double E_rel; + Direction v_target; + while (true) { + // sample target velocity with the constant cross section (cxs) approx. + v_target = sample_cxs_target_velocity(nuc->awr_, E, u, kT); + Direction v_rel = v_neut - v_target; + E_rel = v_rel.dot(v_rel); + if (E_rel < E_up) break; + } -// // get max 0K xs value over range of practical relative energies -// xs_max = std::max(xs_low, & -// maxval(nuc % elastic_0K(i_E_low + 1 : i_E_up)), xs_up) + // perform Doppler broadening rejection correction (dbrc) + double xs_0K = nuc->elastic_xs_0K(E_rel); + double R = xs_0K / xs_max; + if (prn() < R) return v_target; + } -// DBRC_REJECT_LOOP: do -// TARGET_ENERGY_LOOP: do -// // sample target velocity with the constant cross section (cxs) approx. -// sample_cxs_target_velocity(nuc, v_target, E, uvw, kT) -// E_rel = dot_product((v_neut - v_target), (v_neut - v_target)) -// if (E_rel < E_up) exit TARGET_ENERGY_LOOP -// end do TARGET_ENERGY_LOOP + } else if (sampling_method == RES_SCAT_ARES) { + // interpolate xs CDF since we're not exactly at the energy indices + // cdf value at lower bound attainable energy + double m = (nuc->xs_cdf_[i_E_low] - nuc->xs_cdf_[i_E_low - 1]) + / (nuc->energy_0K_[i_E_low + 1] - nuc->energy_0K_[i_E_low]); + double cdf_low = nuc->xs_cdf_[i_E_low - 1] + + m * (E_low - nuc->energy_0K_[i_E_low]); + if (E_low <= nuc->energy_0K_.front()) cdf_low = 0.0; -// // perform Doppler broadening rejection correction (dbrc) -// xs_0K = elastic_xs_0K(E_rel, nuc) -// R = xs_0K / xs_max -// if (prn() < R) exit DBRC_REJECT_LOOP -// end do DBRC_REJECT_LOOP + // cdf value at upper bound attainable energy + m = (nuc->xs_cdf_[i_E_up] - nuc->xs_cdf_[i_E_up - 1]) + / (nuc->energy_0K_[i_E_up + 1] - nuc->energy_0K_[i_E_up]); + double cdf_up = nuc->xs_cdf_[i_E_up - 1] + + m*(E_up - nuc->energy_0K_[i_E_up]); -// } else if (sampling_method == RES_SCAT_ARES) { -// // interpolate xs CDF since we're not exactly at the energy indices -// // cdf value at lower bound attainable energy -// m = (nuc % xs_cdf(i_E_low) - nuc % xs_cdf(i_E_low - 1)) & -// / (nuc % energy_0K(i_E_low + 1) - nuc % energy_0K(i_E_low)) -// cdf_low = nuc % xs_cdf(i_E_low - 1) & -// + m * (E_low - nuc % energy_0K(i_E_low)) -// if (E_low <= nuc % energy_0K(1)) cdf_low = 0.0 + while (true) { + // directly sample Maxwellian + double E_t = -kT * std::log(prn()); -// // cdf value at upper bound attainable energy -// m = (nuc % xs_cdf(i_E_up) - nuc % xs_cdf(i_E_up - 1)) & -// / (nuc % energy_0K(i_E_up + 1) - nuc % energy_0K(i_E_up)) -// cdf_up = nuc % xs_cdf(i_E_up - 1) & -// + m * (E_up - nuc % energy_0K(i_E_up)) + // sample a relative energy using the xs cdf + double cdf_rel = cdf_low + prn()*(cdf_up - cdf_low); + int i_E_rel = lower_bound_index(&nuc->xs_cdf_[i_E_low-1], + &nuc->xs_cdf_[i_E_up+1], cdf_rel); + double E_rel = nuc->energy_0K_[i_E_low + i_E_rel]; + double m = (nuc->xs_cdf_[i_E_low + i_E_rel] + - nuc->xs_cdf_[i_E_low + i_E_rel - 1]) + / (nuc->energy_0K_[i_E_low + i_E_rel + 1] + - nuc->energy_0K_[i_E_low + i_E_rel]); + E_rel += (cdf_rel - nuc->xs_cdf_[i_E_low + i_E_rel - 1]) / m; -// ARES_REJECT_LOOP: do + // perform rejection sampling on cosine between + // neutron and target velocities + double mu = (E_t + nuc->awr_ * (E - E_rel)) / + (2.0 * std::sqrt(nuc->awr_ * E * E_t)); -// // directly sample Maxwellian -// E_t = -kT * std::log(prn()) + if (std::abs(mu) < 1.0) { + // set and accept target velocity + E_t /= nuc->awr_; + return std::sqrt(E_t) * rotate_angle(u, mu, nullptr); + } + } + } + } + } +} -// // sample a relative energy using the xs cdf -// cdf_rel = cdf_low + prn() * (cdf_up - cdf_low) -// i_E_rel = binary_search(nuc % xs_cdf(i_E_low-1:i_E_up), & -// i_E_up - i_E_low + 2, cdf_rel) -// E_rel = nuc % energy_0K(i_E_low + i_E_rel - 1) -// m = (nuc % xs_cdf(i_E_low + i_E_rel - 1) & -// - nuc % xs_cdf(i_E_low + i_E_rel - 2)) & -// / (nuc % energy_0K(i_E_low + i_E_rel) & -// - nuc % energy_0K(i_E_low + i_E_rel - 1)) -// E_rel = E_rel + (cdf_rel - nuc % xs_cdf(i_E_low + i_E_rel - 2)) / m +Direction +sample_cxs_target_velocity(double awr, double E, Direction u, double kT) +{ + double beta_vn = std::sqrt(awr * E / kT); + double alpha = 1.0/(1.0 + std::sqrt(PI)*beta_vn/2.0); -// // perform rejection sampling on cosine between -// // neutron and target velocities -// mu = (E_t + awr * (E - E_rel)) / (2.0 * std::sqrt(awr * E * E_t)) + double beta_vt_sq; + double mu; + while (true) { + // Sample two random numbers + double r1 = prn(); + double r2 = prn(); -// if (abs(mu) < 1.0) { -// // set and accept target velocity -// E_t = E_t / awr -// v_target = std::sqrt(E_t) * rotate_angle(uvw, mu) -// exit ARES_REJECT_LOOP -// } -// end do ARES_REJECT_LOOP -// } -// } -// end select -// } + if (prn() < alpha) { + // With probability alpha, we sample the distribution p(y) = + // y*e^(-y). This can be done with sampling scheme C45 frmo the Monte + // Carlo sampler -// void sample_cxs_target_velocity(int i_nuclide, Direction* v_target, double E, Direction u, -// double kT) -// { -// awr = nuc % awr + beta_vt_sq = -std::log(r1*r2); -// beta_vn = std::sqrt(awr * E / kT) -// alpha = 1.0/(1.0 + std::sqrt(pi)*beta_vn/2.0) + } else { + // With probability 1-alpha, we sample the distribution p(y) = y^2 * + // e^(-y^2). This can be done with sampling scheme C61 from the Monte + // Carlo sampler -// do -// // Sample two random numbers -// r1 = prn() -// r2 = prn() + double c = std::cos(PI/2.0 * prn()); + beta_vt_sq = -std::log(r1) - std::log(r2)*c*c; + } -// if (prn() < alpha) { -// // With probability alpha, we sample the distribution p(y) = -// // y*e^(-y). This can be done with sampling scheme C45 frmo the Monte -// // Carlo sampler + // Determine beta * vt + double beta_vt = std::sqrt(beta_vt_sq); -// beta_vt_sq = -std::log(r1*r2) + // Sample cosine of angle between neutron and target velocity + mu = 2.0*prn() - 1.0; -// } else { -// // With probability 1-alpha, we sample the distribution p(y) = y^2 * -// // e^(-y^2). This can be done with sampling scheme C61 from the Monte -// // Carlo sampler + // Determine rejection probability + double accept_prob = std::sqrt(beta_vn*beta_vn + beta_vt_sq - + 2*beta_vn*beta_vt*mu) / (beta_vn + beta_vt); -// c = std::cos(PI/2.0 * prn()) -// beta_vt_sq = -std::log(r1) - std::log(r2)*c*c -// } + // Perform rejection sampling on vt and mu + if (prn() < accept_prob) break; + } -// // Determine beta * vt -// beta_vt = std::sqrt(beta_vt_sq) + // Determine speed of target nucleus + double vt = std::sqrt(beta_vt_sq*kT/awr); -// // Sample cosine of angle between neutron and target velocity -// mu = 2.0*prn() - 1.0 - -// // Determine rejection probability -// accept_prob = std::sqrt(beta_vn*beta_vn + beta_vt_sq - 2*beta_vn*beta_vt*mu) & -// /(beta_vn + beta_vt) - -// // Perform rejection sampling on vt and mu -// if (prn() < accept_prob) exit -// end do - -// // Determine speed of target nucleus -// vt = std::sqrt(beta_vt_sq*kT/awr) - -// // Determine velocity vector of target nucleus based on neutron's velocity -// // and the sampled angle between them -// v_target = vt * rotate_angle(uvw, mu) -// } + // Determine velocity vector of target nucleus based on neutron's velocity + // and the sampled angle between them + return vt * rotate_angle(u, mu, nullptr); +} void sample_fission_neutron(int i_nuclide, const Reaction* rx, double E_in, Bank* site) { @@ -1073,53 +1092,55 @@ void sample_fission_neutron(int i_nuclide, const Reaction* rx, double E_in, Bank } } -// void inelastic_scatter(int i_nuclide, const Reaction& rx, Particle* p) -// { -// // copy energy of neutron -// E_in = p->E +void inelastic_scatter(const Nuclide* nuc, const Reaction* rx, Particle* p) +{ + // copy energy of neutron + double E_in = p->E; -// // sample outgoing energy and scattering cosine -// rxn % product_sample(1, E_in, E, mu) + // sample outgoing energy and scattering cosine + double E; + double mu; + rx->products_[0].sample(E_in, E, mu); -// // if scattering system is in center-of-mass, transfer cosine of scattering -// // angle and outgoing energy from CM to LAB -// if (rxn % scatter_in_cm) { -// E_cm = E + // if scattering system is in center-of-mass, transfer cosine of scattering + // angle and outgoing energy from CM to LAB + if (rx->scatter_in_cm_) { + double E_cm = E; -// // determine outgoing energy in lab -// A = nuc%awr -// E = E_cm + (E_in + 2.0 * mu * (A+1.0) * std::sqrt(E_in * E_cm)) & -// / ((A+1.0)*(A+1.0)) + // determine outgoing energy in lab + double A = nuc->awr_; + E = E_cm + (E_in + 2.0*mu*(A + 1.0) * std::sqrt(E_in*E_cm)) + / ((A + 1.0)*(A + 1.0)); -// // determine outgoing angle in lab -// mu = mu * std::sqrt(E_cm/E) + 1.0/(A+1.0) * std::sqrt(E_in/E) -// } + // determine outgoing angle in lab + mu = mu*std::sqrt(E_cm/E) + 1.0/(A+1.0) * std::sqrt(E_in/E); + } -// // Because of floating-point roundoff, it may be possible for mu to be -// // outside of the range [-1,1). In these cases, we just set mu to exactly -1 -// // or 1 -// if (abs(mu) > 1.0) mu = sign(1.0,mu) + // Because of floating-point roundoff, it may be possible for mu to be + // outside of the range [-1,1). In these cases, we just set mu to exactly -1 + // or 1 + if (std::abs(mu) > 1.0) mu = std::copysign(1.0, mu); -// // Set outgoing energy and scattering angle -// p->E = E -// p->mu = mu + // Set outgoing energy and scattering angle + p->E = E; + p->mu = mu; -// // change direction of particle -// p->coord(1) % uvw = rotate_angle(p->coord(1) % uvw, mu) + // change direction of particle + rotate_angle_c(p->coord[0].uvw, mu, nullptr); -// // evaluate yield -// yield = rxn % product_yield(1, E_in) -// if (mod(yield, 1.0) == 0.0) { -// // If yield is integral, create exactly that many secondary particles -// do i = 1, nint(yield) - 1 -// particle_create_secondary(p, p->coord(1) % uvw, p->E, & -// NEUTRON, run_CE=true) -// end do -// } else { -// // Otherwise, change weight of particle based on yield -// p->wgt = yield * p->wgt -// } -// } + // evaluate yield + double yield = (*rx->products_[0].yield_)(E_in); + if (std::floor(yield) == yield) { + // If yield is integral, create exactly that many secondary particles + for (int i = 0; i < static_cast(std::round(yield)) - 1; ++i) { + int neutron = static_cast(ParticleType::neutron); + p->create_secondary(p->coord[0].uvw, p->E, neutron, true); + } + } else { + // Otherwise, change weight of particle based on yield + p->wgt *= yield; + } +} void sample_secondary_photons(Particle* p, int i_nuclide) { diff --git a/src/position.cpp b/src/position.cpp index 18ec33e803..ecf5be631f 100644 --- a/src/position.cpp +++ b/src/position.cpp @@ -60,4 +60,22 @@ Position::operator*=(double v) return *this; } +Position& +Position::operator/=(Position other) +{ + x /= other.x; + y /= other.y; + z /= other.z; + return *this; +} + +Position& +Position::operator/=(double v) +{ + x /= v; + y /= v; + z /= v; + return *this; +} + } // namespace openmc diff --git a/src/settings.cpp b/src/settings.cpp index 36657c6c62..303a35ddd9 100644 --- a/src/settings.cpp +++ b/src/settings.cpp @@ -806,18 +806,6 @@ void read_settings_xml() //============================================================================== extern "C" { - bool res_scat_nuclides_empty() { - return settings::res_scat_nuclides.empty(); - } - - int res_scat_nuclides_size() { - return settings::res_scat_nuclides.size(); - } - - bool res_scat_nuclides_cmp(int i, const char* name) { - return settings::res_scat_nuclides[i - 1] == name; - } - const char* path_cross_sections_c() { return settings::path_cross_sections.c_str(); }