From 68c135856513630e2fedf4873db8e72a8397c715 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Wed, 11 Jul 2018 13:21:05 -0500 Subject: [PATCH] Use Reaction type hierarchy on C++ side --- src/endf.cpp | 5 + src/endf.h | 1 + src/nuclide_header.F90 | 146 ++++++++--------- src/physics.F90 | 56 ++++--- src/reaction.cpp | 115 +++++++++++++- src/reaction.h | 23 +++ src/reaction_header.F90 | 293 ++++++++++++++++++++++++++++------- src/reaction_product.h | 2 +- src/secondary_uncorrelated.h | 4 +- src/tallies/tally.F90 | 49 +++--- 10 files changed, 502 insertions(+), 192 deletions(-) diff --git a/src/endf.cpp b/src/endf.cpp index ed23cc336..c0c84acfa 100644 --- a/src/endf.cpp +++ b/src/endf.cpp @@ -28,6 +28,11 @@ Interpolation int2interp(int i) } } +bool is_fission(int mt) +{ + return mt == 18 || mt == 19 || mt == 20 || mt == 21 || mt == 38; +} + //============================================================================== // Polynomial implementation //============================================================================== diff --git a/src/endf.h b/src/endf.h index e3fce4dfb..189b919b1 100644 --- a/src/endf.h +++ b/src/endf.h @@ -9,6 +9,7 @@ namespace openmc { Interpolation int2interp(int i); +bool is_fission(int MT); class Function1D { public: diff --git a/src/nuclide_header.F90 b/src/nuclide_header.F90 index 931df199e..319bf93ab 100644 --- a/src/nuclide_header.F90 +++ b/src/nuclide_header.F90 @@ -642,29 +642,33 @@ contains if (rx % MT >= N_2N0 .and. rx % MT <= N_2NC .and. find(MTs, N_2N) /= -1) cycle do t = 1, n_temperature - j = rx % xs(t) % threshold - n = size(rx % xs(t) % value) + j = rx % xs_threshold(t) + n = rx % xs_size(t) ! Add contribution to total cross section - this % xs(t) % value(XS_TOTAL,j:j+n-1) = this % xs(t) % & - value(XS_TOTAL,j:j+n-1) + rx % xs(t) % value + do k = j, j + n - 1 + this % xs(t) % value(XS_TOTAL,k) = this % xs(t) % & + value(XS_TOTAL,k) + rx % xs(t, k - j + 1) + end do ! Calculate photon production cross section - do k = 1, size(rx % products) - if (rx % products(k) % particle == PHOTON) then + do k = 1, rx % products_size() + if (rx % product_particle(k) == PHOTON) then do l = 1, n this % xs(t) % value(XS_PHOTON_PROD,l+j-1) = & this % xs(t) % value(XS_PHOTON_PROD,l+j-1) + & - rx % xs(t) % value(l) * rx % products(k) % & - yield % evaluate(this % grid(t) % energy(l+j-1)) + rx % xs(t, l) * rx % product_yield(k, & + this % grid(t) % energy(l+j-1)) end do end if end do ! Add contribution to absorption cross section if (is_disappearance(rx % MT)) then - this % xs(t) % value(XS_ABSORPTION,j:j+n-1) = this % xs(t) % & - value(XS_ABSORPTION,j:j+n-1) + rx % xs(t) % value + do k = j, j + n - 1 + this % xs(t) % value(XS_ABSORPTION,k) = this % xs(t) % & + value(XS_ABSORPTION,k) + rx % xs(t, k - j + 1) + end do end if ! Information about fission reactions @@ -680,39 +684,20 @@ contains ! Add contribution to fission cross section if (is_fission(rx % MT)) then this % fissionable = .true. - this % xs(t) % value(XS_FISSION,j:j+n-1) = this % xs(t) % & - value(XS_FISSION,j:j+n-1) + rx % xs(t) % value + do k = j, j + n - 1 + this % xs(t) % value(XS_FISSION,k) = this % xs(t) % & + value(XS_FISSION,k) + rx % xs(t, k - j + 1) - ! Also need to add fission cross sections to absorption - this % xs(t) % value(XS_ABSORPTION,j:j+n-1) = this % xs(t) % & - value(XS_ABSORPTION,j:j+n-1) + rx % xs(t) % value + ! Also need to add fission cross sections to absorption + this % xs(t) % value(XS_ABSORPTION,k) = this % xs(t) % & + value(XS_ABSORPTION,k) + rx % xs(t, k - j + 1) + end do ! Keep track of this reaction for easy searching later if (t == 1) then i_fission = i_fission + 1 this % index_fission(i_fission) = i this % n_fission = this % n_fission + 1 - - ! <<<<<<<<<<<<<<<<<<<<<<<<<<<< REMOVE THIS <<<<<<<<<<<<<<<<<<<<<<<<< - ! Before the secondary distribution refactor, when the angle/energy - ! distribution was uncorrelated, no angle was actually sampled. With - ! the refactor, an angle is always sampled for an uncorrelated - ! distribution even when no angle distribution exists in the ACE file - ! (isotropic is assumed). To preserve the RNG stream, we explicitly - ! mark fission reactions so that we avoid the angle sampling. - do k = 1, size(rx % products) - if (rx % products(k) % particle == NEUTRON) then - do m = 1, size(rx % products(k) % distribution) - associate (aedist => rx % products(k) % distribution(m) % obj) - select type (aedist) - type is (UncorrelatedAngleEnergy) - aedist % fission = .true. - end select - end associate - end do - end if - end do - ! <<<<<<<<<<<<<<<<<<<<<<<<<<<< REMOVE THIS <<<<<<<<<<<<<<<<<<<<<<<<< end if end if ! fission end do ! temperature @@ -721,12 +706,13 @@ contains ! Determine number of delayed neutron precursors if (this % fissionable) then - do i = 1, size(this % reactions(this % index_fission(1)) % products) - if (this % reactions(this % index_fission(1)) % products(i) % & - emission_mode == EMISSION_DELAYED) then - this % n_precursor = this % n_precursor + 1 - end if - end do + associate (rx => this % reactions(this % index_fission(1))) + do i = 1, rx % products_size() + if (rx % product_emission_mode(i) == EMISSION_DELAYED) then + this % n_precursor = this % n_precursor + 1 + end if + end do + end associate end if ! Calculate nu-fission cross section @@ -761,36 +747,30 @@ contains select case (emission_mode) case (EMISSION_PROMPT) - associate (product => this % reactions(this % index_fission(1)) % products(1)) - nu = product % yield % evaluate(E) + associate (rx => this % reactions(this % index_fission(1))) + nu = rx % product_yield(1, E) end associate case (EMISSION_DELAYED) if (this % n_precursor > 0) then - if (present(group) .and. group < & - size(this % reactions(this % index_fission(1)) % products)) then - ! If delayed group specified, determine yield immediately - associate(p => this % reactions(this % index_fission(1)) % products(1 + group)) - nu = p % yield % evaluate(E) - end associate + associate(rx => this % reactions(this % index_fission(1))) + if (present(group) .and. group < rx % products_size()) then + ! If delayed group specified, determine yield immediately + nu = rx % product_yield(1 + group, E) + else + nu = ZERO - else - nu = ZERO + do i = 2, rx % products_size() + ! Skip any non-neutron products + if (rx % product_particle(i) /= NEUTRON) exit - associate (rx => this % reactions(this % index_fission(1))) - do i = 2, size(rx % products) - associate (product => rx % products(i)) - ! Skip any non-neutron products - if (product % particle /= NEUTRON) exit - - ! Evaluate yield - if (product % emission_mode == EMISSION_DELAYED) then - nu = nu + product % yield % evaluate(E) - end if - end associate + ! Evaluate yield + if (rx % product_emission_mode(i) == EMISSION_DELAYED) then + nu = nu + rx % product_yield(i, E) + end if end do - end associate - end if + end if + end associate else nu = ZERO end if @@ -799,8 +779,8 @@ contains if (allocated(this % total_nu)) then nu = this % total_nu % evaluate(E) else - associate (product => this % reactions(this % index_fission(1)) % products(1)) - nu = product % yield % evaluate(E) + associate (rx => this % reactions(this % index_fission(1))) + nu = rx % product_yield(1, E) end associate end if end select @@ -868,6 +848,7 @@ contains 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_t, sig_a, sig_f ! Intermediate multipole variables @@ -1009,10 +990,11 @@ contains ! need to specifically check its threshold index i_rxn = this % reaction_index(DEPLETION_RX(1)) if (i_rxn > 0) then - associate (xs => this % reactions(i_rxn) % xs(i_temp)) + associate (rx => this % reactions(i_rxn)) + threshold = rx % xs_threshold(i_temp) micro_xs % reaction(1) = (ONE - f) * & - xs % value(i_grid - xs % threshold + 1) + & - f * xs % value(i_grid - xs % threshold + 2) + rx % xs(i_temp, i_grid - threshold + 1) + & + f * rx % xs(i_temp, i_grid - threshold + 2) end associate end if @@ -1022,11 +1004,12 @@ contains ! reaction xs appropriately i_rxn = this % reaction_index(DEPLETION_RX(j)) if (i_rxn > 0) then - associate (xs => this % reactions(i_rxn) % xs(i_temp)) - if (i_grid >= xs % threshold) then + associate (rx => this % reactions(i_rxn)) + threshold = rx % xs_threshold(i_temp) + if (i_grid >= threshold) then micro_xs % reaction(j) = (ONE - f) * & - xs % value(i_grid - xs % threshold + 1) + & - f * xs % value(i_grid - xs % threshold + 2) + 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 @@ -1091,8 +1074,9 @@ contains f = micro_xs % interp_factor if (i_temp > 0) then - associate (xs => this % reactions(1) % xs(i_temp) % value) - micro_xs % elastic = (ONE - f) * xs(i_grid) + f * xs(i_grid + 1) + associate (rx => this % reactions(1)) + micro_xs % elastic = (ONE - f) * rx % xs(i_temp, i_grid) + & + f * rx % xs(i_temp, i_grid + 1) end associate else ! For multipole, elastic is total - absorption @@ -1459,6 +1443,7 @@ contains integer :: i_energy ! index for energy integer :: i_low ! band index at lower bounding energy integer :: i_up ! band index at upper bounding energy + integer :: threshold ! threshold energy index real(8) :: f ! interpolation factor real(8) :: r ! pseudo-random number real(8) :: elastic ! elastic cross section @@ -1551,10 +1536,11 @@ contains f = micro_xs % interp_factor ! Determine inelastic scattering cross section - associate (xs => this % reactions(this % urr_inelastic) % xs(i_temp)) - if (i_energy >= xs % threshold) then - inelastic = (ONE - f) * xs % value(i_energy - xs % threshold + 1) + & - f * xs % value(i_energy - xs % threshold + 2) + associate (rx => this % reactions(this % urr_inelastic)) + threshold = rx % xs_threshold(i_temp) + if (i_energy >= threshold) then + inelastic = (ONE - f) * rx % xs(i_temp, i_energy - threshold + 1) + & + f * rx % xs(i_temp, i_energy - threshold + 2) end if end associate end if diff --git a/src/physics.F90 b/src/physics.F90 index 0eed83418..dab6bb7f9 100644 --- a/src/physics.F90 +++ b/src/physics.F90 @@ -506,6 +506,7 @@ contains integer :: i integer :: i_grid integer :: i_temp + integer :: threshold real(8) :: f real(8) :: prob real(8) :: cutoff @@ -545,13 +546,14 @@ contains FISSION_REACTION_LOOP: do i = 1, nuc % n_fission i_reaction = nuc % index_fission(i) - associate (xs => nuc % reactions(i_reaction) % xs(i_temp)) + associate (rx => nuc % reactions(i_reaction)) ! if energy is below threshold for this reaction, skip it - if (i_grid < xs % threshold) cycle + threshold = rx % xs_threshold(i_temp) + if (i_grid < threshold) cycle ! add to cumulative probability - prob = prob + ((ONE - f) * xs % value(i_grid - xs % threshold + 1) & - + f*(xs % value(i_grid - xs % threshold + 2))) + prob = prob + ((ONE - f) * rx % xs(i_temp, i_grid - threshold + 1) & + + f*(rx % xs(i_temp, i_grid - threshold + 2))) end associate ! Create fission bank sites if fission occurs @@ -593,17 +595,17 @@ contains ! Loop through each reaction type REACTION_LOOP: do i_reaction = 1, size(nuc % reactions) associate (rx => nuc % reactions(i_reaction)) - threshold = rx % xs(i_temp) % threshold + threshold = rx % xs_threshold(i_temp) ! if energy is below threshold for this reaction, skip it if (i_grid < threshold) cycle - do i_product = 1, size(rx % products) - if (rx % products(i_product) % particle == PHOTON) then + do i_product = 1, rx % products_size() + if (rx % product_particle(i_product) == PHOTON) then ! add to cumulative probability - yield = rx % products(i_product) % yield % evaluate(E) - prob = prob + ((ONE - f) * rx % xs(i_temp) % value(i_grid - threshold + 1) & - + f*(rx % xs(i_temp) % value(i_grid - threshold + 2))) * yield + yield = rx % product_yield(i_product, E) + prob = prob + ((ONE - f) * rx % xs(i_temp, i_grid - threshold + 1) & + + f*(rx % xs(i_temp, i_grid - threshold + 2))) * yield if (prob > cutoff) return last_valid_reaction = i_reaction @@ -672,6 +674,7 @@ contains integer :: j integer :: i_temp integer :: i_grid + integer :: threshold real(8) :: f real(8) :: prob real(8) :: cutoff @@ -750,14 +753,14 @@ contains &// trim(nuc % name)) end if - associate (rx => nuc % reactions(i), & - xs => nuc % reactions(i) % xs(i_temp)) + associate (rx => nuc % reactions(i)) ! if energy is below threshold for this reaction, skip it - if (i_grid < xs % threshold) cycle + threshold = rx % xs_threshold(i_temp) + if (i_grid < threshold) cycle ! add to cumulative probability - prob = prob + ((ONE - f)*xs % value(i_grid - xs % threshold + 1) & - + f*(xs % value(i_grid - xs % threshold + 2))) + 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 @@ -839,14 +842,7 @@ contains vel = sqrt(dot_product(v_n, v_n)) ! Sample scattering angle - select type (dist => rxn % products(1) % distribution(1) % obj) - type is (UncorrelatedAngleEnergy) - if (allocated(dist % angle % energy)) then - mu_cm = dist % angle % sample(E) - else - mu_cm = TWO*prn() - ONE - end if - end select + mu_cm = rxn % sample_elastic_mu(E) ! Determine direction cosines in CM uvw_cm = v_n/vel @@ -1581,7 +1577,7 @@ contains do group = 1, nuc % n_precursor ! determine delayed neutron precursor yield for group j - yield = rxn % products(1 + group) % yield % evaluate(E_in) + yield = rxn % product_yield(1 + group, E_in) ! Check if this group is sampled prob = prob + yield @@ -1600,7 +1596,7 @@ contains do ! sample from energy/angle distribution -- note that mu has already been ! sampled above and doesn't need to be resampled - call rxn % products(1 + group) % sample(E_in, site % E, mu) + call rxn % product_sample(1 + group, E_in, site % E, mu) ! resample if energy is greater than maximum neutron energy if (site % E < energy_max(NEUTRON)) exit @@ -1624,7 +1620,7 @@ contains ! sample from prompt neutron energy distribution n_sample = 0 do - call rxn % products(1) % sample(E_in, site % E, mu) + call rxn % product_sample(1, E_in, site % E, mu) ! resample if energy is greater than maximum neutron energy if (site % E < energy_max(NEUTRON)) exit @@ -1663,7 +1659,7 @@ contains E_in = p % E ! sample outgoing energy and scattering cosine - call rxn % products(1) % sample(E_in, E, mu) + 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 @@ -1692,7 +1688,7 @@ contains p % coord(1) % uvw = rotate_angle(p % coord(1) % uvw, mu) ! evaluate yield - yield = rxn % products(1) % yield % evaluate(E_in) + 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 @@ -1740,8 +1736,8 @@ contains call sample_photon_product(i_nuclide, p % E, i_reaction, i_product) ! Sample the outgoing energy and angle - call nuclides(i_nuclide) % reactions(i_reaction) % products(i_product) & - % sample(p % E, E, mu) + call nuclides(i_nuclide) % reactions(i_reaction) % & + product_sample(i_product, p % E, E, mu) ! Sample the new direction uvw = rotate_angle(p % coord(1) % uvw, mu) diff --git a/src/reaction.cpp b/src/reaction.cpp index 7e6a3ae5d..5e89155ae 100644 --- a/src/reaction.cpp +++ b/src/reaction.cpp @@ -4,6 +4,9 @@ #include // for move #include "hdf5_interface.h" +#include "endf.h" +#include "random_lcg.h" +#include "secondary_uncorrelated.h" namespace openmc { @@ -43,6 +46,116 @@ Reaction::Reaction(hid_t group, const std::vector& temperatures) close_group(pgroup); } } + + // <<<<<<<<<<<<<<<<<<<<<<<<<<<< REMOVE THIS <<<<<<<<<<<<<<<<<<<<<<<<< + // Before the secondary distribution refactor, when the angle/energy + // distribution was uncorrelated, no angle was actually sampled. With + // the refactor, an angle is always sampled for an uncorrelated + // distribution even when no angle distribution exists in the ACE file + // (isotropic is assumed). To preserve the RNG stream, we explicitly + // mark fission reactions so that we avoid the angle sampling. + if (is_fission(mt_)) { + for (auto& p : products_) { + if (p.particle_ == ParticleType::neutron) { + for (auto& d : p.distribution_) { + auto d_ = dynamic_cast(d.get()); + if (d_) d_->fission_ = true; + } + } + } + } + // <<<<<<<<<<<<<<<<<<<<<<<<<<<< REMOVE THIS <<<<<<<<<<<<<<<<<<<<<<<<< } -} \ No newline at end of file +//============================================================================== +// Fortran compatibility functions +//============================================================================== + +Reaction* reaction_from_hdf5(hid_t group, int* temperatures, int n) +{ + std::vector temps {temperatures, temperatures + n}; + return new Reaction{group, temps}; +} + +void reaction_delete(Reaction* rx) { delete rx; } + +int reaction_mt(Reaction* rx) { return rx->mt_; } + +double reaction_q_value(Reaction* rx) { return rx->q_value_; } + +bool reaction_scatter_in_cm(Reaction* rx) { return rx->scatter_in_cm_; } + +double reaction_product_decay_rate(Reaction* rx, int product) +{ + return rx->products_[product - 1].decay_rate_; +} + +int reaction_product_emission_mode(Reaction* rx, int product) +{ + switch (rx->products_[product - 1].emission_mode_) { + case ReactionProduct::EmissionMode::prompt: + return 1; + case ReactionProduct::EmissionMode::delayed: + return 2; + case ReactionProduct::EmissionMode::total: + return 3; + } +} + +int reaction_product_particle(Reaction* rx, int product) +{ + switch (rx->products_[product - 1].particle_) { + case ParticleType::neutron: + return 1; + case ParticleType::photon: + return 2; + case ParticleType::electron: + return 3; + case ParticleType::positron: + return 4; + } +} + +void reaction_product_sample(Reaction* rx, int product, double E_in, double* E_out, double* mu) +{ + rx->products_[product - 1].sample(E_in, *E_out, *mu); +} + +double reaction_product_yield(Reaction* rx, int product, double E) +{ + return (*rx->products_[product - 1].yield_)(E); +} + +int reaction_products_size(Reaction* rx) { return rx->products_.size(); } + +double reaction_xs(Reaction* rx, int temperature, int energy) +{ + return rx->xs_[temperature - 1].value[energy - 1]; +} + +double reaction_sample_elastic_mu(Reaction* rx, double E) +{ + // Get elastic scattering distribution + auto& d = rx->products_[0].distribution_[0]; + + // Check if it is an uncorrelated angle-energy distribution + auto d_ = dynamic_cast(d.get()); + if (d_) { + return d_->sampleMu(E); + } else { + return 2.0*prn() - 1.0; + } + +} + +int reaction_xs_size(Reaction* rx, int temperature) +{ + return rx->xs_[temperature - 1].value.size(); +} + +int reaction_xs_threshold(Reaction* rx, int temperature) +{ + return rx->xs_[temperature - 1].threshold; +} + +} diff --git a/src/reaction.h b/src/reaction.h index 0040bd094..170554edf 100644 --- a/src/reaction.h +++ b/src/reaction.h @@ -24,6 +24,29 @@ public: std::vector products_; }; +//============================================================================== +// Fortran compatibility functions +//============================================================================== + +extern "C" { + Reaction* reaction_from_hdf5(hid_t group, int* temperatures, int n); + void reaction_delete(Reaction* rx); + int reaction_mt(Reaction* rx); + double reaction_q_value(Reaction* rx); + bool reaction_scatter_in_cm(Reaction* rx); + double reaction_product_decay_rate(Reaction* rx, int product); + int reaction_product_emission_mode(Reaction* rx, int product); + int reaction_product_particle(Reaction* rx, int product); + void reaction_product_sample(Reaction* rx, int product, double E_in, + double* E_out, double* mu); + int reaction_products_size(Reaction* rx); + double reaction_product_yield(Reaction* rx, int product, double E); + double reaction_sample_elastic_mu(Reaction* rx, double E); + double reaction_xs(Reaction* xs, int temperature, int energy); + int reaction_xs_size(Reaction* xs, int temperature); + int reaction_xs_threshold(Reaction* xs, int temperature); +} + } #endif // OPENMC_REACTION_H diff --git a/src/reaction_header.F90 b/src/reaction_header.F90 index 48973336b..915bf4833 100644 --- a/src/reaction_header.F90 +++ b/src/reaction_header.F90 @@ -1,83 +1,268 @@ module reaction_header + use, intrinsic :: ISO_C_BINDING + use constants, only: MAX_WORD_LEN use hdf5_interface - use product_header, only: ReactionProduct use stl_vector, only: VectorInt use string, only: to_str, starts_with implicit none + private !=============================================================================== ! REACTION contains the cross-section and secondary energy and angle ! distributions for a single reaction in a continuous-energy ACE-format table !=============================================================================== - type TemperatureXS - integer :: threshold ! Energy grid index of threshold - real(8), allocatable :: value(:) ! Cross section values - end type TemperatureXS - - type Reaction - integer :: MT ! ENDF MT value - real(8) :: Q_value ! Reaction Q value - logical :: scatter_in_cm ! scattering system in center-of-mass? - type(TemperatureXS), allocatable :: xs(:) - type(ReactionProduct), allocatable :: products(:) + type, public :: Reaction + type(C_PTR) :: ptr + integer(C_INT) :: MT ! ENDF MT value + real(C_DOUBLE) :: Q_value ! Reaction Q value + logical(C_BOOL) :: scatter_in_cm ! scattering system in center-of-mass? contains - procedure :: from_hdf5 => reaction_from_hdf5 + procedure :: from_hdf5 + procedure :: mt_ + procedure :: q_value_ + procedure :: scatter_in_cm_ + procedure :: product_decay_rate + procedure :: product_emission_mode + procedure :: product_particle + procedure :: product_sample + procedure :: product_yield + procedure :: products_size + procedure :: sample_elastic_mu + procedure :: xs + procedure :: xs_size + procedure :: xs_threshold end type Reaction + interface + function reaction_from_hdf5(group, temperatures, n) result(ptr) bind(C) + import C_PTR, HID_T, C_INT + integer(HID_T), value :: group + integer(C_INT), intent(in) :: temperatures + integer(C_INT), value :: n + type(C_PTR) :: ptr + end function + + function reaction_mt(ptr) result(mt) bind(C) + import C_PTR, C_INT + type(C_PTR), value :: ptr + integer(C_INT) :: mt + end function + + function reaction_q_value(ptr) result(q_value) bind(C) + import C_PTR, C_DOUBLE + type(C_PTR), value :: ptr + real(C_DOUBLE) :: q_value + end function + + function reaction_scatter_in_cm(ptr) result(b) bind(C) + import C_PTR, C_BOOL + type(C_PTR), value :: ptr + logical(C_BOOL) :: b + end function + + pure function reaction_product_decay_rate(ptr, product) result(rate) bind(C) + import C_PTR, C_INT, C_DOUBLE + type(C_PTR), value :: ptr + integer(C_INT), value :: product + real(C_DOUBLE) :: rate + end function + + pure function reaction_product_emission_mode(ptr, product) result(m) bind(C) + import C_PTR, C_INT + type(C_PTR), value :: ptr + integer(C_INT), value :: product + integer(C_INT) :: m + end function + + pure function reaction_product_particle(ptr, product) result(particle) bind(C) + import C_PTR, C_INT + type(C_PTR), value :: ptr + integer(C_INT), value :: product + integer(C_INT) :: particle + end function + + subroutine reaction_product_sample(ptr, product, E_in, E_out, mu) bind(C) + import C_PTR, C_INT, C_DOUBLE + type(C_PTR), value :: ptr + integer(C_INT), value :: product + real(C_DOUBLE), value :: E_in + real(C_DOUBLE), intent(out) :: E_out + real(C_DOUBLE), intent(out) :: mu + end subroutine + + pure function reaction_product_yield(ptr, product, E) result(val) bind(C) + import C_PTR, C_INT, C_DOUBLE + type(C_PTR), value :: ptr + integer(C_INT), value :: product + real(C_DOUBLE), value :: E + real(C_DOUBLE) :: val + end function + + pure function reaction_products_size(ptr) result(sz) bind(C) + import C_PTR, C_INT + type(C_PTR), value :: ptr + integer(C_INT) :: sz + end function + + function reaction_sample_elastic_mu(ptr, E) result(mu) bind(C) + import C_PTR, C_INT, C_DOUBLE + type(C_PTR), value :: ptr + real(C_DOUBLE), value :: E + real(C_DOUBLE) :: mu + end function + + function reaction_xs(ptr, temperature, energy) result(xs) bind(C) + import C_PTR, C_INT, C_DOUBLE + type(C_PTR), value :: ptr + integer(C_INT), value :: temperature + integer(C_INT), value :: energy + real(C_DOUBLE) :: xs + end function + + function reaction_xs_size(ptr, temperature) result(sz) bind(C) + import C_PTR, C_INT + type(C_PTR), value :: ptr + integer(C_INT), value :: temperature + integer(C_INT) :: sz + end function + + function reaction_xs_threshold(ptr, temperature) result(threshold) bind(C) + import C_PTR, C_INT + type(C_PTR), value :: ptr + integer(C_INT), value :: temperature + integer(C_INT) :: threshold + end function + end interface + contains - subroutine reaction_from_hdf5(this, group_id, temperatures) + subroutine from_hdf5(this, group_id, temperatures) class(Reaction), intent(inout) :: this integer(HID_T), intent(in) :: group_id type(VectorInt), intent(in) :: temperatures - integer :: i - integer :: cm - integer :: n_product - integer(HID_T) :: pgroup - integer(HID_T) :: xs, temp_group - integer(HSIZE_T) :: dims(1) - integer(HSIZE_T) :: j - character(MAX_WORD_LEN) :: temp_str ! temperature dataset name, e.g. '294K' - character(MAX_WORD_LEN), allocatable :: grp_names(:) + integer(C_INT) :: dummy + integer(C_INT) :: n - call read_attribute(this % Q_value, group_id, 'Q_value') - call read_attribute(this % MT, group_id, 'mt') - call read_attribute(cm, group_id, 'center_of_mass') - this % scatter_in_cm = (cm == 1) + n = temperatures % size() + if (n > 0) then + this % ptr = reaction_from_hdf5(group_id, temperatures % data(1), n) + else + ! In this case, temperatures % data(1) doesn't exist, so we just pass a + ! dummy value + this % ptr = reaction_from_hdf5(group_id, dummy, n) + end if + this % MT = reaction_mt(this % ptr) + this % Q_value = reaction_q_value(this % ptr) + this % scatter_in_cm = reaction_scatter_in_cm(this % ptr) + end subroutine from_hdf5 - ! Read cross section and threshold_idx data - allocate(this % xs(temperatures % size())) - do i = 1, temperatures % size() - temp_str = trim(to_str(temperatures % data(i))) // "K" - temp_group = open_group(group_id, temp_str) - xs = open_dataset(temp_group, 'xs') - call read_attribute(this % xs(i) % threshold, xs, 'threshold_idx') - call get_shape(xs, dims) - allocate(this % xs(i) % value(dims(1))) - call read_dataset(this % xs(i) % value, xs) - call close_dataset(xs) - call close_group(temp_group) - end do + function mt_(this) result(mt) + class(Reaction), intent(in) :: this + integer(C_INT) :: MT - ! Determine number of products - n_product = 0 - call get_groups(group_id, grp_names) - do j = 1, size(grp_names) - if (starts_with(grp_names(j), "product_")) n_product = n_product + 1 - end do + mt = reaction_mt(this % ptr) + end function - ! Read products - allocate(this % products(n_product)) - do i = 1, n_product - pgroup = open_group(group_id, 'product_' // trim(to_str(i - 1))) - call this % products(i) % from_hdf5(pgroup) - call close_group(pgroup) - end do - end subroutine reaction_from_hdf5 + function q_value_(this) result(q_value) + class(Reaction), intent(in) :: this + real(C_DOUBLE) :: q_value + + q_value = reaction_q_value(this % ptr) + end function + + function scatter_in_cm_(this) result(cm) + class (Reaction), intent(in) :: this + logical(C_BOOL) :: cm + + cm = reaction_scatter_in_cm(this % ptr) + end function + + pure function product_decay_rate(this, product) result(rate) + class(Reaction), intent(in) :: this + integer(C_INT), intent(in) :: product + real(C_DOUBLE) :: rate + + rate = reaction_product_decay_rate(this % ptr, product) + end function + + pure function product_emission_mode(this, product) result(m) + class(Reaction), intent(in) :: this + integer(C_INT), intent(in) :: product + integer(C_INT) :: m + + m = reaction_product_emission_mode(this % ptr, product) + end function + + pure function product_particle(this, product) result(p) + class(Reaction), intent(in) :: this + integer(C_INT), intent(in) :: product + integer(C_INT) :: p + + p = reaction_product_particle(this % ptr, product) + end function + + subroutine product_sample(this, product, E_in, E_out, mu) + class(Reaction), intent(in) :: this + integer(C_INT), intent(in) :: product + real(C_DOUBLE), intent(in) :: E_in + real(C_DOUBLE), intent(out) :: E_out + real(C_DOUBLE), intent(out) :: mu + + call reaction_product_sample(this % ptr, product, E_in, E_out, mu) + end subroutine + + pure function product_yield(this, product, E) result(val) + class(Reaction), intent(in) :: this + integer(C_INT), intent(in) :: product + real(C_DOUBLE), intent(in) :: E + real(C_DOUBLE) :: val + + val = reaction_product_yield(this % ptr, product, E) + end function + + pure function products_size(this) result(sz) + class(Reaction), intent(in) :: this + integer(C_INT) :: sz + + sz = reaction_products_size(this % ptr) + end function + + function sample_elastic_mu(this, E) result(mu) + class(Reaction), intent(in) :: this + real(C_DOUBLE), intent(in) :: E + real(C_DOUBLE) :: mu + + mu = reaction_sample_elastic_mu(this % ptr, E) + end function + + function xs(this, temperature, energy) result(val) + class(Reaction), intent(in) :: this + integer(C_INT), intent(in) :: temperature + integer(C_INT), intent(in) :: energy + real(C_DOUBLE) :: val + + val = reaction_xs(this % ptr, temperature, energy) + end function + + function xs_size(this, temperature) result(sz) + class(Reaction), intent(in) :: this + integer(C_INT) :: temperature + integer(C_INT) :: sz + + sz = reaction_xs_size(this % ptr, temperature) + end function + + function xs_threshold(this, temperature) result(val) + class(Reaction), intent(in) :: this + integer(C_INT), intent(in) :: temperature + integer(C_INT) :: val + + val = reaction_xs_threshold(this % ptr, temperature) + end function end module reaction_header diff --git a/src/reaction_product.h b/src/reaction_product.h index c72657f00..2ad55866b 100644 --- a/src/reaction_product.h +++ b/src/reaction_product.h @@ -15,7 +15,7 @@ class ReactionProduct { public: explicit ReactionProduct(hid_t group); void sample(double E_in, double& E_out, double& mu) const; -private: + enum class EmissionMode { prompt, // Prompt emission of secondary particle total, // Delayed emission of secondary particle diff --git a/src/secondary_uncorrelated.h b/src/secondary_uncorrelated.h index 73c62c14a..3205fb8fd 100644 --- a/src/secondary_uncorrelated.h +++ b/src/secondary_uncorrelated.h @@ -15,8 +15,10 @@ class UncorrelatedAngleEnergy : public AngleEnergy { public: explicit UncorrelatedAngleEnergy(hid_t group); void sample(double E_in, double& E_out, double& mu) const; -private: + double sampleMu(double E) const { return angle_.sample(E); } + bool fission_ {false}; +private: AngleDistribution angle_; std::unique_ptr energy_; }; diff --git a/src/tallies/tally.F90 b/src/tallies/tally.F90 index 8ddd42bce..2c660ce87 100644 --- a/src/tallies/tally.F90 +++ b/src/tallies/tally.F90 @@ -94,6 +94,7 @@ contains integer :: k ! loop index for bank sites integer :: d_bin ! delayed group bin index integer :: dg_filter ! index of delayed group filter + integer :: threshold ! threshold energy index real(8) :: yield ! delayed neutron yield real(8) :: atom_density_ ! atom/b-cm real(8) :: f ! interpolation factor @@ -227,7 +228,7 @@ contains ! Get yield and apply to score associate (rxn => nuclides(p % event_nuclide) % reactions(m)) score = p % last_wgt * flux & - * rxn % products(1) % yield % evaluate(E) + * rxn % product_yield(1, E) end associate end if @@ -633,7 +634,7 @@ contains score = p % absorb_wgt * yield * & micro_xs(p % event_nuclide) % fission & / micro_xs(p % event_nuclide) % absorption & - * rxn % products(1 + d) % decay_rate * flux + * rxn % product_decay_rate(1 + d) * flux end associate ! Tally to bin @@ -657,9 +658,8 @@ contains ! rxn % products array to be exceeded. Hence, we use the size ! of this array and not the MAX_DELAYED_GROUPS constant for ! this loop. - do d = 1, size(rxn % products) - 2 - - score = score + rxn % products(1 + d) % decay_rate * & + do d = 1, rxn % products_size() - 2 + score = score + rxn % product_decay_rate(1 + d) * & p % absorb_wgt & * micro_xs(p % event_nuclide) % fission & * nuclides(p % event_nuclide) % & @@ -699,7 +699,7 @@ contains ! determine score based on bank site weight and keff. score = score + keff * fission_bank(n_bank - p % n_bank + k) & - % wgt * rxn % products(1 + g) % decay_rate * flux + % wgt * rxn % product_decay_rate(1 + g) * flux end associate ! if the delayed group filter is present, tally to corresponding @@ -755,7 +755,7 @@ contains ! Compute the score and tally to bin score = micro_xs(i_nuclide) % fission * yield * flux * & - atom_density * rxn % products(1 + d) % decay_rate + atom_density * rxn % product_decay_rate(1 + d) end associate ! Tally to bin @@ -778,11 +778,10 @@ contains ! groups since this could cause the range of the rxn % products ! array to be exceeded. Hence, we use the size of this array ! and not the MAX_DELAYED_GROUPS constant for this loop. - do d = 1, size(rxn % products) - 2 - + do d = 1, rxn % products_size() - 2 score = score + micro_xs(i_nuclide) % fission * flux * & nuclides(i_nuclide) % nu(E, EMISSION_DELAYED) * & - atom_density * rxn % products(1 + d) % decay_rate + atom_density * rxn % product_decay_rate(1 + d) end do end associate end if @@ -824,7 +823,7 @@ contains ! Compute the score score = micro_xs(i_nuc) % fission * yield * flux * & atom_density_ & - * rxn % products(1 + d) % decay_rate + * rxn % product_decay_rate(1 + d) end associate ! Tally to bin @@ -860,13 +859,13 @@ contains ! rxn % products array to be exceeded. Hence, we use the ! size of this array and not the MAX_DELAYED_GROUPS ! constant for this loop. - do d = 1, size(rxn % products) - 2 + do d = 1, rxn % products_size() - 2 ! Accumulate the contribution from each nuclide score = score + micro_xs(i_nuc) % fission & * nuclides(i_nuc) % nu(E, EMISSION_DELAYED) & * atom_density_ * flux & - * rxn % products(1 + d) % decay_rate + * rxn % product_decay_rate(1 + d) end do end associate end if @@ -1132,12 +1131,12 @@ contains i_energy = micro_xs(i_nuclide) % index_grid f = micro_xs(i_nuclide) % interp_factor - associate (xs => nuclides(i_nuclide) % reactions(m) & - % xs(i_temp)) - if (i_energy >= xs % threshold) then - score = ((ONE - f) * xs % value(i_energy - & - xs % threshold + 1) + f * xs % value(i_energy - & - xs % threshold + 2)) * atom_density * flux + associate (rx => nuclides(i_nuclide) % reactions(m)) + threshold = rx % xs_threshold(i_temp) + if (i_energy >= threshold) then + score = ((ONE - f) * rx % xs(i_temp, i_energy - & + threshold + 1) + f * rx % xs(i_temp, i_energy - & + threshold + 2)) * atom_density * flux end if end associate else @@ -1165,12 +1164,12 @@ contains i_energy = micro_xs(i_nuc) % index_grid f = micro_xs(i_nuc) % interp_factor - associate (xs => nuclides(i_nuc) % reactions(m) & - % xs(i_temp)) - if (i_energy >= xs % threshold) then - score = score + ((ONE - f) * xs % value(i_energy - & - xs % threshold + 1) + f * xs % value(i_energy - & - xs % threshold + 2)) * atom_density_ * flux + associate (rx => nuclides(i_nuc) % reactions(m)) + threshold = rx % xs_threshold(i_temp) + if (i_energy >= threshold) then + score = score + ((ONE - f) * rx % xs(i_temp, i_energy - & + threshold + 1) + f * rx % xs(i_temp, i_energy - & + threshold + 2)) * atom_density_ * flux end if end associate else