diff --git a/include/openmc/physics.h b/include/openmc/physics.h index 58bc7e990..1b5fb9023 100644 --- a/include/openmc/physics.h +++ b/include/openmc/physics.h @@ -58,7 +58,7 @@ void scatter(Particle*, int i_nuclide, int i_nuc_mat); void elastic_scatter(int i_nuclide, const Reaction* rx, double kT, double* E, double* uvw, double* mu_lab, double* wgt); -extern "C" void sab_scatter(int i_nuclide, int i_sab, double* E, +void sab_scatter(int i_nuclide, int i_sab, double* E, double* uvw, double* mu); //! samples the target velocity. The constant cross section free gas model is diff --git a/include/openmc/thermal.h b/include/openmc/thermal.h index 265831e10..66596b76d 100644 --- a/include/openmc/thermal.h +++ b/include/openmc/thermal.h @@ -2,6 +2,7 @@ #define OPENMC_THERMAL_H #include +#include #include #include @@ -27,6 +28,16 @@ constexpr int SAB_SECONDARY_CONT {2}; // Continuous, linear-linear interpolati constexpr int SAB_ELASTIC_INCOHERENT {3}; // Incoherent elastic scattering constexpr int SAB_ELASTIC_COHERENT {4}; // Coherent elastic scattering (Bragg edges) +//============================================================================== +// Global variables +//============================================================================== + +class ThermalScattering; + +namespace data { +extern std::vector> thermal_scatt; +} + //============================================================================== //! Secondary angle-energy data for thermal neutron scattering at a single //! temperature @@ -37,7 +48,7 @@ public: ThermalData(hid_t group, int secondary_mode); // Sample an outgoing energy and angle - void sample(const NuclideMicroXS* micro_xs, double E_in, + void sample(const NuclideMicroXS& micro_xs, double E_in, double* E_out, double* mu); private: //! Secondary energy/angle distributions for inelastic thermal scattering @@ -112,7 +123,7 @@ public: bool has_nuclide(const char* name) const; // Sample an outgoing energy and angle - void sample(const NuclideMicroXS* micro_xs, double E_in, + void sample(const NuclideMicroXS& micro_xs, double E_in, double* E_out, double* mu); double threshold() const { return data_[0].threshold_inelastic_; } @@ -133,12 +144,8 @@ public: extern "C" { ThermalScattering* sab_from_hdf5(hid_t group, const double* temperature, int n, int method, double tolerance, const double* minmax); - void sab_calculate_xs(ThermalScattering* data, double E, double sqrtkT, - int* i_temp, double* elastic, double* inelastic); void sab_free(ThermalScattering* data); bool sab_has_nuclide(ThermalScattering* data, const char* name); - void sab_sample(ThermalScattering* data, const NuclideMicroXS* micro_xs, - double E_in, double* E_out, double* mu); double sab_threshold(ThermalScattering* data); } diff --git a/src/input_xml.F90 b/src/input_xml.F90 index 41bac0f88..f0796c3f6 100644 --- a/src/input_xml.F90 +++ b/src/input_xml.F90 @@ -2157,10 +2157,6 @@ contains end if ! Set up logarithmic grid for nuclides - do i = 1, size(nuclides) - call nuclides(i) % init_grid(energy_min(NEUTRON), & - energy_max(NEUTRON), n_log_bins) - end do log_spacing = log(energy_max(NEUTRON)/energy_min(NEUTRON)) / n_log_bins do i = 1, size(materials) diff --git a/src/material_header.F90 b/src/material_header.F90 index e041f52f4..13914cc59 100644 --- a/src/material_header.F90 +++ b/src/material_header.F90 @@ -446,7 +446,7 @@ contains .or. i_sab /= micro_xs(i_nuclide) % index_sab & .or. sab_frac /= micro_xs(i_nuclide) % sab_frac) then call nuclides(i_nuclide) % calculate_xs(i_sab, p % E, i_grid, & - p % sqrtkT, sab_frac, micro_xs(i_nuclide)) + p % sqrtkT, sab_frac) end if ! ====================================================================== diff --git a/src/nuclide.cpp b/src/nuclide.cpp index 188847677..011364588 100644 --- a/src/nuclide.cpp +++ b/src/nuclide.cpp @@ -10,6 +10,7 @@ #include "openmc/settings.h" #include "openmc/simulation.h" #include "openmc/string_utils.h" +#include "openmc/thermal.h" #include "xtensor/xbuilder.hpp" #include "xtensor/xview.hpp" @@ -683,7 +684,13 @@ void Nuclide::calculate_xs(int i_sab, double E, int i_log_union, // If the particle is in the unresolved resonance range and there are // probability tables, we need to determine cross sections from the table - this->calculate_urr_xs(i_temp, E); + if (settings::urr_ptables_on && urr_present_ && !use_mp) { + int n = urr_data_[i_temp].n_energy_; + if ((E > urr_data_[i_temp].energy_(0)) && + (E < urr_data_[i_temp].energy_(n-1))) { + this->calculate_urr_xs(i_temp, E); + } + } micro_xs.last_E = E; micro_xs.last_sqrtkT = sqrtkT; @@ -700,7 +707,7 @@ void Nuclide::calculate_sab_xs(int i_sab, double E, double sqrtkT, double sab_fr int i_temp; double elastic; double inelastic; - //data::sab_tables[i_sab]->calculate_xs(E, sqrtkT, &i_temp, &elastic, &inelastic); + data::thermal_scatt[i_sab]->calculate_xs(E, sqrtkT, &i_temp, &elastic, &inelastic); // Store the S(a,b) cross sections. micro.thermal = sab_frac * (elastic + inelastic); @@ -869,6 +876,12 @@ extern "C" Reaction* nuclide_reaction(Nuclide* nuc, int i_rx) return nuc->reactions_[i_rx-1].get(); } +extern "C" void nuclide_calculate_xs_c(Nuclide* nuc, int i_sab, double E, + int i_log_union, double sqrtkT, double sab_frac) +{ + nuc->calculate_xs(i_sab, E, i_log_union, sqrtkT, sab_frac); +} + extern "C" void nuclides_clear() { data::nuclides.clear(); } @@ -881,17 +894,4 @@ void set_micro_xs() } } -extern "C" void -nuclide_calculate_urr_xs(bool use_mp, int i_nuclide, int i_temp, double E) -{ - Nuclide* nuc = data::nuclides[i_nuclide - 1].get(); - if (settings::urr_ptables_on && (nuc->urr_present_ && !use_mp)) { - if ((E > nuc->urr_data_[i_temp - 1].energy_(0)) && - (E < nuc->urr_data_[i_temp - 1].energy_( - nuc->urr_data_[i_temp - 1].n_energy_ - 1))) { - nuc->calculate_urr_xs(i_temp - 1, E); - } - } -} - } // namespace openmc diff --git a/src/nuclide_header.F90 b/src/nuclide_header.F90 index fe35f7293..04c543522 100644 --- a/src/nuclide_header.F90 +++ b/src/nuclide_header.F90 @@ -32,7 +32,6 @@ module nuclide_header !=============================================================================== type EnergyGrid - integer, allocatable :: grid_index(:) ! log grid mapping indices real(8), allocatable :: energy(:) ! energy values corresponding to xs end type EnergyGrid @@ -95,7 +94,6 @@ module nuclide_header contains procedure :: clear => nuclide_clear procedure :: from_hdf5 => nuclide_from_hdf5 - procedure :: init_grid => nuclide_init_grid procedure :: nu => nuclide_nu procedure, private :: create_derived => nuclide_create_derived procedure :: calculate_xs => nuclide_calculate_xs @@ -192,14 +190,6 @@ module nuclide_header character(kind=C_CHAR), intent(in) :: name(*) type(C_PTR) :: path end function - - subroutine nuclide_calculate_urr_xs(use_mp, i_nuclide, i_temp, E) bind(C) - import C_BOOL, C_INT, C_DOUBLE - logical(C_BOOL), value, intent(in) :: use_mp - integer(C_INT), value, intent(in) :: i_nuclide - integer(C_INT), value, intent(in) :: i_temp - real(C_DOUBLE), value, intent(in) :: E - end subroutine nuclide_calculate_urr_xs end interface contains @@ -683,51 +673,13 @@ contains end function nuclide_nu - subroutine nuclide_init_grid(this, E_min, E_max, M) - class(Nuclide), intent(inout) :: this - real(8), intent(in) :: E_min ! Minimum energy in eV - real(8), intent(in) :: E_max ! Maximum energy in eV - integer, intent(in) :: M ! Number of equally log-spaced bins - - integer :: i, j, k ! Loop indices - integer :: t ! temperature index - real(8) :: spacing - real(8), allocatable :: umesh(:) ! Equally log-spaced energy grid - - ! Determine equal-logarithmic energy spacing - spacing = log(E_max/E_min)/M - - ! Create equally log-spaced energy grid - allocate(umesh(0:M)) - umesh(:) = [(i*spacing, i=0, M)] - - do t = 1, size(this % grid) - ! Allocate logarithmic mapping for nuclide - allocate(this % grid(t) % grid_index(0:M)) - - ! Determine corresponding indices in nuclide grid to energies on - ! equal-logarithmic grid - j = 1 - do k = 0, M - do while (log(this % grid(t) % energy(j + 1)/E_min) <= umesh(k)) - ! Ensure that for isotopes where maxval(this % energy) << E_max - ! that there are no out-of-bounds issues. - if (j + 1 == size(this % grid(t) % energy)) exit - j = j + 1 - end do - this % grid(t) % grid_index(k) = j - end do - end do - - end subroutine nuclide_init_grid - !=============================================================================== ! NUCLIDE_CALCULATE_XS determines microscopic cross sections for the nuclide ! at the energy of the given particle !=============================================================================== subroutine nuclide_calculate_xs(this, i_sab, E, i_log_union, sqrtkT, & - sab_frac, micro_xs) + sab_frac) class(Nuclide), intent(in) :: this ! Nuclide object integer, intent(in) :: i_sab ! index into sab_tables array real(8), intent(in) :: E ! energy @@ -735,214 +687,19 @@ contains ! material union energy grid real(8), intent(in) :: sqrtkT ! square root of kT, material dependent real(8), intent(in) :: sab_frac ! fraction of atoms affected by S(a,b) - type(NuclideMicroXS), intent(inout) :: micro_xs ! Cross section cache - logical(C_BOOL) :: use_mp ! true if XS can be calculated with windowed multipole - integer :: i_temp ! index for temperature - integer :: i_grid ! index on nuclide energy grid - integer :: i_low ! lower logarithmic mapping index - integer :: i_high ! upper logarithmic mapping index - integer :: i_rxn ! reaction index - integer :: j ! index in DEPLETION_RX - integer :: threshold ! threshold energy index - real(8) :: f ! interp factor on nuclide energy grid - real(8) :: kT ! temperature in eV - real(8) :: sig_s, sig_a, sig_f ! Intermediate multipole variables - - ! Initialize cached cross sections to zero - micro_xs % elastic = CACHE_INVALID - micro_xs % thermal = ZERO - micro_xs % thermal_elastic = ZERO - - ! Check to see if there is multipole data present at this energy - use_mp = .false. - if (this % mp_present) then - if (E >= this % multipole % E_min .and. & - E <= this % multipole % E_max) then - use_mp = .true. - end if - end if - - ! Evaluate multipole or interpolate - if (use_mp) then - ! Call multipole kernel - call multipole_eval(this % multipole, E, sqrtkT, sig_s, sig_a, sig_f) - - micro_xs % total = sig_s + sig_a - micro_xs % elastic = sig_s - micro_xs % absorption = sig_a - micro_xs % fission = sig_f - - if (this % fissionable) then - micro_xs % nu_fission = sig_f * this % nu(E, EMISSION_TOTAL) - else - micro_xs % nu_fission = ZERO - end if - - if (need_depletion_rx) then - ! Initialize all reaction cross sections to zero - micro_xs % reaction(:) = ZERO - - ! Only non-zero reaction is (n,gamma) - micro_xs % reaction(1) = sig_a - sig_f - end if - - ! Ensure these values are set - ! Note, the only time either is used is in one of 4 places: - ! 1. physics.F90 - scatter - For inelastic scatter. - ! 2. physics.F90 - sample_fission - For partial fissions. - ! 3. tally.F90 - score_general - For tallying on MTxxx reactions. - ! 4. nuclide.h - calculate_urr_xs - For unresolved purposes. - ! It is worth noting that none of these occur in the resolved - ! resonance range, so the value here does not matter. index_temp is - ! set to -1 to force a segfault in case a developer messes up and tries - ! to use it with multipole. - micro_xs % index_temp = -1 - micro_xs % index_grid = 0 - micro_xs % interp_factor = ZERO - - else - ! Find the appropriate temperature index. - kT = sqrtkT**2 - select case (temperature_method) - case (TEMPERATURE_NEAREST) - i_temp = minloc(abs(this % kTs - kT), dim=1) - - case (TEMPERATURE_INTERPOLATION) - ! Find temperatures that bound the actual temperature - do i_temp = 1, size(this % kTs) - 1 - if (this % kTs(i_temp) <= kT .and. kT < this % kTs(i_temp + 1)) exit - end do - - ! Randomly sample between temperature i and i+1 - f = (kT - this % kTs(i_temp)) / & - (this % kTs(i_temp + 1) - this % kTs(i_temp)) - if (f > prn()) i_temp = i_temp + 1 - end select - - associate (grid => this % grid(i_temp), xs => this % xs(i_temp)) - ! Determine the energy grid index using a logarithmic mapping to - ! reduce the energy range over which a binary search needs to be - ! performed - - if (E < grid % energy(1)) then - i_grid = 1 - elseif (E > grid % energy(size(grid % energy))) then - i_grid = size(grid % energy) - 1 - else - ! Determine bounding indices based on which equal log-spaced - ! interval the energy is in - i_low = grid % grid_index(i_log_union) - i_high = grid % grid_index(i_log_union + 1) + 1 - - ! Perform binary search over reduced range - i_grid = binary_search(grid % energy(i_low:i_high), & - i_high - i_low + 1, E) + i_low - 1 - end if - - ! check for rare case where two energy points are the same - if (grid % energy(i_grid) == grid % energy(i_grid + 1)) & - i_grid = i_grid + 1 - - ! calculate interpolation factor - f = (E - grid % energy(i_grid)) / & - (grid % energy(i_grid + 1) - grid % energy(i_grid)) - - micro_xs % index_temp = i_temp - micro_xs % index_grid = i_grid - micro_xs % interp_factor = f - - ! Calculate microscopic nuclide total cross section - micro_xs % total = (ONE - f) * xs % value(XS_TOTAL,i_grid) & - + f * xs % value(XS_TOTAL,i_grid + 1) - - ! Calculate microscopic nuclide absorption cross section - micro_xs % absorption = (ONE - f) * xs % value(XS_ABSORPTION, & - i_grid) + f * xs % value(XS_ABSORPTION,i_grid + 1) - - if (this % fissionable) then - ! Calculate microscopic nuclide total cross section - micro_xs % fission = (ONE - f) * xs % value(XS_FISSION,i_grid) & - + f * xs % value(XS_FISSION,i_grid + 1) - - ! Calculate microscopic nuclide nu-fission cross section - micro_xs % nu_fission = (ONE - f) * xs % value(XS_NU_FISSION, & - i_grid) + f * xs % value(XS_NU_FISSION,i_grid + 1) - else - micro_xs % fission = ZERO - micro_xs % nu_fission = ZERO - end if - - ! Calculate microscopic nuclide photon production cross section - micro_xs % photon_prod = (ONE - f) * xs % value(XS_PHOTON_PROD,i_grid) & - + f * xs % value(XS_PHOTON_PROD,i_grid + 1) - end associate - - ! Depletion-related reactions - if (need_depletion_rx) then - ! Initialize all reaction cross sections to zero - micro_xs % reaction(:) = ZERO - - ! Physics says that (n,gamma) is not a threshold reaction, so we don't - ! need to specifically check its threshold index - i_rxn = this % reaction_index(DEPLETION_RX(1)) - if (i_rxn > 0) then - associate (rx => this % reactions(i_rxn)) - threshold = rx % xs_threshold(i_temp) - micro_xs % reaction(1) = (ONE - f) * & - rx % xs(i_temp, i_grid - threshold + 1) + & - f * rx % xs(i_temp, i_grid - threshold + 2) - end associate - end if - - ! Loop over remaining depletion reactions - do j = 2, 6 - ! If reaction is present and energy is greater than threshold, set the - ! reaction xs appropriately - i_rxn = this % reaction_index(DEPLETION_RX(j)) - if (i_rxn > 0) then - associate (rx => this % reactions(i_rxn)) - threshold = rx % xs_threshold(i_temp) - if (i_grid >= threshold) then - micro_xs % reaction(j) = (ONE - f) * & - rx % xs(i_temp, i_grid - threshold + 1) + & - f * rx % xs(i_temp, i_grid - threshold + 2) - elseif (j >= 4) then - ! One can show that the the threshold for (n,(x+1)n) is always - ! higher than the threshold for (n,xn). Thus, if we are below - ! the threshold for, e.g., (n,2n), there is no reason to check - ! the threshold for (n,3n) and (n,4n). - exit - end if - end associate - end if - end do - end if - end if - - ! Initialize sab treatment to false - micro_xs % index_sab = NONE - micro_xs % sab_frac = ZERO - - ! Initialize URR probability table treatment to false - micro_xs % use_ptable = .false. - - ! If there is S(a,b) data for this nuclide, we need to set the sab_scatter - ! and sab_elastic cross sections and correct the total and elastic cross - ! sections. - - if (i_sab > 0) then - call calculate_sab_xs(this, i_sab, E, sqrtkT, sab_frac, micro_xs) - end if - - - ! If the particle is in the unresolved resonance range and there are - ! probability tables, we need to determine cross sections from the table - call nuclide_calculate_urr_xs(use_mp, this % i_nuclide, i_temp, E) - - micro_xs % last_E = E - micro_xs % last_sqrtkT = sqrtkT + interface + subroutine nuclide_calculate_xs_c(i_sab, E, i_log_union, sqrtkT, sab_frac) bind(C) + import C_INT, C_DOUBLE + integer(C_INT), value :: i_sab + real(C_DOUBLE), value :: E + integer(C_INT), value :: i_log_union + real(C_DOUBLE), value :: sqrtkT + real(C_DOUBLE), value :: sab_frac + end subroutine + end interface + call nuclide_calculate_xs_c(i_sab, E, i_log_union, sqrtkT, sab_frac) end subroutine nuclide_calculate_xs !=============================================================================== @@ -973,49 +730,6 @@ contains end if end subroutine nuclide_calculate_elastic_xs -!=============================================================================== -! CALCULATE_SAB_XS determines the elastic and inelastic scattering -! cross-sections in the thermal energy range. These cross sections replace a -! fraction of whatever data were taken from the normal Nuclide table. -!=============================================================================== - - subroutine calculate_sab_xs(this, i_sab, E, sqrtkT, sab_frac, micro_xs) - class(Nuclide), intent(in) :: this ! Nuclide object - integer, intent(in) :: i_sab ! index into sab_tables array - real(8), intent(in) :: E ! energy - real(8), intent(in) :: sqrtkT ! temperature - real(8), intent(in) :: sab_frac ! fraction of atoms affected by S(a,b) - type(NuclideMicroXS), intent(inout) :: micro_xs ! Cross section cache - - integer(C_INT) :: i_temp ! temperature index - real(C_DOUBLE) :: inelastic ! S(a,b) inelastic cross section - real(C_DOUBLE) :: elastic ! S(a,b) elastic cross section - - ! Set flag that S(a,b) treatment should be used for scattering - micro_xs % index_sab = i_sab - - ! Calculate the S(a,b) cross section - call sab_tables(i_sab) % calculate_xs(E, sqrtkT, i_temp, elastic, inelastic) - - ! Store the S(a,b) cross sections. - micro_xs % thermal = sab_frac * (elastic + inelastic) - micro_xs % thermal_elastic = sab_frac * elastic - - ! Calculate free atom elastic cross section - call this % calculate_elastic_xs(micro_xs) - - ! Correct total and elastic cross sections - micro_xs % total = micro_xs % total + micro_xs % thermal - & - sab_frac * micro_xs % elastic - micro_xs % elastic = micro_xs % thermal + (ONE - sab_frac) * & - micro_xs % elastic - - ! Save temperature index and thermal fraction - micro_xs % index_temp_sab = i_temp - micro_xs % sab_frac = sab_frac - - end subroutine calculate_sab_xs - !=============================================================================== ! MULTIPOLE_EVAL evaluates the windowed multipole equations for cross ! sections in the resolved resonance regions @@ -1333,10 +1047,6 @@ contains ! Add entry to nuclide dictionary call nuclide_dict % set(to_lower(name_), n) n_nuclides = n - - ! Initialize nuclide grid - call nuclides(n) % init_grid(energy_min(NEUTRON), & - energy_max(NEUTRON), n_log_bins) else err = E_DATA call set_errmsg("Nuclide '" // trim(name_) // "' is not present & diff --git a/src/physics.F90 b/src/physics.F90 index 0576e2a88..0812592ea 100644 --- a/src/physics.F90 +++ b/src/physics.F90 @@ -249,28 +249,4 @@ contains end function sample_element -!=============================================================================== -! 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) 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 - - ! Sample from C++ side - ptr = C_LOC(micro_xs(i_nuclide + 1)) - call sab_tables(i_sab) % sample(ptr, E, E_out, mu) - - ! Set energy to outgoing, change direction of particle - E = E_out - uvw = rotate_angle(uvw, mu) - end subroutine sab_scatter - end module physics diff --git a/src/physics.cpp b/src/physics.cpp index 27d0e7b81..e3af8bdef 100644 --- a/src/physics.cpp +++ b/src/physics.cpp @@ -782,16 +782,20 @@ void elastic_scatter(int i_nuclide, const Reaction* rx, double kT, double* E, 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) -// { -// // Sample from C++ side -// ptr = C_LOC(micro_xs(i_nuclide)) -// sab_tables(i_sab) % sample(ptr, E, E_out, mu) +void sab_scatter(int i_nuclide, int i_sab, double* E, double* uvw, double* mu) +{ + // Determine temperature index + const auto& micro {simulation::micro_xs[i_nuclide]}; + int i_temp = micro.index_temp_sab; -// // Set energy to outgoing, change direction of particle -// E = E_out -// uvw = rotate_angle(uvw, mu) -// } + // Sample energy and angle + double E_out; + data::thermal_scatt[i_sab]->data_[i_temp].sample(micro, *E, &E_out, mu); + + // Set energy to outgoing, change direction of particle + *E = E_out; + rotate_angle_c(uvw, *mu, nullptr); +} Direction sample_target_velocity(const Nuclide* nuc, double E, Direction u, Direction v_neut, double xs_eff, double kT, double* wgt) diff --git a/src/sab_header.F90 b/src/sab_header.F90 index b4696a65d..061da5d1c 100644 --- a/src/sab_header.F90 +++ b/src/sab_header.F90 @@ -21,10 +21,8 @@ module sab_header type(C_PTR) :: ptr contains procedure :: from_hdf5 - procedure :: calculate_xs procedure :: free procedure :: has_nuclide - procedure :: sample procedure :: threshold end type SAlphaBeta @@ -46,17 +44,6 @@ module sab_header type(C_PTR) :: ptr end function - subroutine sab_calculate_xs(ptr, E, sqrtkT, i_temp, elastic, & - inelastic) bind(C) - import C_PTR, C_DOUBLE, C_INT - type(C_PTR), value :: ptr - real(C_DOUBLE), value :: E - real(C_DOUBLE), value :: sqrtkT - integer(C_INT), intent(out) :: i_temp - real(C_DOUBLE), intent(out) :: elastic - real(C_DOUBLE), intent(out) :: inelastic - end subroutine - subroutine sab_free(ptr) bind(C) import C_PTR type(C_PTR), value :: ptr @@ -69,15 +56,6 @@ module sab_header logical(C_BOOL) :: val end function - subroutine sab_sample(ptr, micro_xs, E_in, E_out, mu) bind(C) - import C_PTR, C_INT, C_DOUBLE - type(C_PTR), value :: ptr - type(C_PTR), value :: micro_xs - real(C_DOUBLE), value :: E_in - real(C_DOUBLE), intent(out) :: E_out - real(C_DOUBLE), intent(out) :: mu - end subroutine - function sab_threshold(ptr) result(threshold) bind(C) import C_PTR, C_double type(C_PTR), value :: ptr @@ -109,22 +87,6 @@ contains end if end subroutine from_hdf5 -!=============================================================================== -! SAB_CALCULATE_XS determines the elastic and inelastic scattering -! cross-sections in the thermal energy range. -!=============================================================================== - - subroutine calculate_xs(this, E, sqrtkT, i_temp, elastic, inelastic) - class(SAlphaBeta), intent(in) :: this ! S(a,b) object - real(C_DOUBLE), intent(in) :: E ! energy - real(C_DOUBLE), intent(in) :: sqrtkT ! temperature - integer, intent(out) :: i_temp ! index in the S(a,b)'s temperature - real(C_DOUBLE), intent(out) :: elastic ! thermal elastic cross section - real(C_DOUBLE), intent(out) :: inelastic ! thermal inelastic cross section - - call sab_calculate_xs(this % ptr, E, sqrtkT, i_temp, elastic, inelastic) - end subroutine - subroutine free(this) class(SAlphaBeta), intent(inout) :: this call sab_free(this % ptr) @@ -138,16 +100,6 @@ contains val = sab_has_nuclide(this % ptr, to_c_string(name)) end function - subroutine sample(this, micro_xs, E_in, E_out, mu) - class(SAlphaBeta), intent(in) :: this - type(C_PTR), value :: micro_xs - real(C_DOUBLE), value :: E_in - real(C_DOUBLE), intent(out) :: E_out - real(C_DOUBLE), intent(out) :: mu - - call sab_sample(this % ptr, micro_xs, E_in, E_out, mu) - end subroutine - function threshold(this) class(SAlphaBeta), intent(in) :: this real(C_DOUBLE) :: threshold diff --git a/src/thermal.cpp b/src/thermal.cpp index 4168d5320..48780e4eb 100644 --- a/src/thermal.cpp +++ b/src/thermal.cpp @@ -20,6 +20,14 @@ namespace openmc { +//============================================================================== +// Global variables +//============================================================================== + +namespace data { +std::vector> thermal_scatt; +} + //============================================================================== // ThermalScattering implementation //============================================================================== @@ -359,11 +367,11 @@ ThermalData::ThermalData(hid_t group, int secondary_mode) } void -ThermalData::sample(const NuclideMicroXS* micro_xs, double E, +ThermalData::sample(const NuclideMicroXS& micro_xs, double E, double* E_out, double* mu) { // Determine whether inelastic or elastic scattering will occur - if (prn() < micro_xs->thermal_elastic / micro_xs->thermal) { + if (prn() < micro_xs.thermal_elastic / micro_xs.thermal) { // elastic scattering // Get index and interpolation factor for elastic grid @@ -590,18 +598,10 @@ sab_from_hdf5(hid_t group, const double* temperature, int n, std::vector T {temperature, temperature + n}; // Create new object and return it - return new ThermalScattering{group, T, method, tolerance, minmax}; -} + data::thermal_scatt.push_back(std::make_unique( + group, T, method, tolerance, minmax)); -void sab_calculate_xs(ThermalScattering* data, double E, double sqrtkT, - int* i_temp, double* elastic, double* inelastic) -{ - // Calculate cross section - int t; - data->calculate_xs(E, sqrtkT, &t, elastic, inelastic); - - // Fortran needs index plus one - *i_temp = t + 1; + return data::thermal_scatt.back().get(); } void sab_free(ThermalScattering* data) { delete data; } @@ -611,13 +611,6 @@ bool sab_has_nuclide(ThermalScattering* data, const char* name) return data->has_nuclide(name); } -void sab_sample(ThermalScattering* data, const NuclideMicroXS* micro_xs, - double E_in, double* E_out, double* mu) -{ - int i_temp = micro_xs->index_temp_sab; - data->data_[i_temp - 1].sample(micro_xs, E_in, E_out, mu); -} - double sab_threshold(ThermalScattering* data) { return data->threshold(); } } // namespace openmc