em/distribution/SBEnergyDistHelper.hh

0001 //----------------------------------*-C++-*----------------------------------//
0002 // Copyright 2021-2024 UT-Battelle, LLC, and other Celeritas developers.
0003 // See the top-level COPYRIGHT file for details.
0004 // SPDX-License-Identifier: (Apache-2.0 OR MIT)
0005 //---------------------------------------------------------------------------//
0006 //! \file celeritas/em/distribution/SBEnergyDistHelper.hh
0007 //---------------------------------------------------------------------------//
0008 #pragma once
0009
0010 #include <cmath>
0011
0012 #include "corecel/grid/TwodGridCalculator.hh"
0013 #include "corecel/grid/TwodSubgridCalculator.hh"
0014 #include "corecel/math/Algorithms.hh"
0015 #include "celeritas/Quantities.hh"
0016 #include "celeritas/em/data/SeltzerBergerData.hh"
0017 #include "celeritas/random/distribution/ReciprocalDistribution.hh"
0018
0019 namespace celeritas
0020 {
0021 //---------------------------------------------------------------------------//
0022 /*!
0023  * Help sample exiting photon energy from Bremsstrahlung.
0024  *
0025  * This class simply preprocesses the input data needed for the
0026  * SBEnergyDistribution, which is templated on a dynamic cross section
0027  * correction factor.
0028  *
0029  * The cross section units are immaterial since the cross section merely acts
0030  * as a shape function for rejection: the sampled energy's cross section is
0031  * always divided by the maximium cross section.
0032  */
0033 class SBEnergyDistHelper
0034 {
0035   public:
0036     //!@{
0037     //! \name Type aliases
0038     using SBDXsec = NativeCRef<SeltzerBergerTableData>;
0039     using Xs = Quantity<SBElementTableData::XsUnits>;
0040     using Energy = units::MevEnergy;
0041     using EnergySq = Quantity<UnitProduct<units::Mev, units::Mev>>;
0042     //!@}
0043
0044   public:
0045     // Construct from data
0046     inline CELER_FUNCTION SBEnergyDistHelper(SBDXsec const& differential_xs,
0047                                              Energy inc_energy,
0048                                              ElementId element,
0049                                              EnergySq density_correction,
0050                                              Energy min_gamma_energy);
0051
0052     // Sample scaled energy (analytic component of exiting distribution)
0053     template<class Engine>
0054     inline CELER_FUNCTION Energy sample_exit_energy(Engine& rng) const;
0055
0056     // Calculate tabulated cross section for a given energy
0057     inline CELER_FUNCTION Xs calc_xs(Energy energy) const;
0058
0059     //! Maximum cross section calculated for rejection
0060     CELER_FUNCTION Xs max_xs() const { return max_xs_; }
0061
0062   private:
0063     //// IMPLEMENTATION TYPES ////
0064
0065     using SBTables = NativeCRef<SeltzerBergerTableData>;
0066     using ReciprocalSampler = ReciprocalDistribution<real_type>;
0067
0068     //// IMPLEMENTATION DATA ////
0069
0070     TwodSubgridCalculator const calc_xs_;
0071     Xs const max_xs_;
0072
0073     real_type const inv_inc_energy_;
0074     real_type const dens_corr_;
0075     ReciprocalSampler const sample_exit_esq_;
0076
0077     //// CONSTRUCTION HELPER FUNCTIONS ////
0078
0079     inline CELER_FUNCTION TwodSubgridCalculator make_xs_calc(
0080         SBTables const&, real_type inc_energy, ElementId element) const;
0081
0082     inline CELER_FUNCTION Xs calc_max_xs(SBTables const& xs_params,
0083                                          ElementId element) const;
0084
0085     inline CELER_FUNCTION ReciprocalSampler
0086     make_esq_sampler(real_type inc_energy, real_type min_gamma_energy) const;
0087 };
0088
0089 //---------------------------------------------------------------------------//
0090 // INLINE DEFINITIONS
0091 //---------------------------------------------------------------------------//
0092 /*!
0093  * Construct from incident particle and energy.
0094  *
0095  * The incident energy *must* be within the bounds of the SB table data, so the
0096  * Model's applicability must be consistent with the table data.
0097  */
0098 CELER_FUNCTION
0099 SBEnergyDistHelper::SBEnergyDistHelper(SBDXsec const& differential_xs,
0100                                        Energy inc_energy,
0101                                        ElementId element,
0102                                        EnergySq density_correction,
0103                                        Energy min_gamma_energy)
0104     : calc_xs_{this->make_xs_calc(differential_xs, inc_energy.value(), element)}
0105     , max_xs_{this->calc_max_xs(differential_xs, element)}
0106     , inv_inc_energy_(1 / inc_energy.value())
0107     , dens_corr_(density_correction.value())
0108     , sample_exit_esq_{
0109           this->make_esq_sampler(inc_energy.value(), min_gamma_energy.value())}
0110 {
0111     CELER_EXPECT(inc_energy > min_gamma_energy);
0112 }
0113
0114 //---------------------------------------------------------------------------//
0115 /*!
0116  * Sample an exit energy on a scaled and adjusted reciprocal distribution.
0117  */
0118 template<class Engine>
0119 CELER_FUNCTION auto
0120 SBEnergyDistHelper::sample_exit_energy(Engine& rng) const -> Energy
0121 {
0122     // Sample scaled energy and subtract correction factor
0123     real_type esq = sample_exit_esq_(rng) - dens_corr_;
0124     CELER_ASSERT(esq >= 0);
0125     return Energy{std::sqrt(esq)};
0126 }
0127
0128 //---------------------------------------------------------------------------//
0129 /*!
0130  * Calculate tabulated cross section for a given energy.
0131  */
0132 CELER_FUNCTION auto SBEnergyDistHelper::calc_xs(Energy e) const -> Xs
0133 {
0134     CELER_EXPECT(e > zero_quantity());
0135     // Interpolate the differential cross setion at the given exit energy
0136     return Xs{calc_xs_(e.value() * inv_inc_energy_)};
0137 }
0138
0139 //---------------------------------------------------------------------------//
0140 /*!
0141  * Construct the differential cross section calculator for exit energy.
0142  */
0143 CELER_FUNCTION TwodSubgridCalculator SBEnergyDistHelper::make_xs_calc(
0144     SBTables const& xs_params, real_type inc_energy, ElementId element) const
0145 {
0146     CELER_EXPECT(element < xs_params.elements.size());
0147     CELER_EXPECT(inc_energy > 0);
0148
0149     TwodGridData const& grid = xs_params.elements[element].grid;
0150     CELER_ASSERT(inc_energy >= std::exp(xs_params.reals[grid.x.front()])
0151                  && inc_energy < std::exp(xs_params.reals[grid.x.back()]));
0152
0153     static_assert(
0154         std::is_same<Energy::unit_type, units::Mev>::value
0155             && std::is_same<SBElementTableData::EnergyUnits, units::LogMev>::value,
0156         "Inconsistent energy units");
0157     return TwodGridCalculator(grid, xs_params.reals)(std::log(inc_energy));
0158 }
0159
0160 //---------------------------------------------------------------------------//
0161 /*!
0162  * Calculate a bounding maximum of the differential cross section.
0163  *
0164  * This interpolates the maximum cross section for the given incident energy
0165  * by using the pre-calculated cross section maxima. The interpolated value is
0166  * typically exactly the
0167  * maximum (since the two \em y points are usually adjacent, and therefore the
0168  * linear interpolation between them is exact) but at worst (e.g. for the
0169  * double-peaked function of brems at lower energies) an upper bound which can
0170  * be proven by the triangle inequality.
0171  *
0172  * \note This is called during construction, so \c calc_xs_ must be initialized
0173  * before whatever calls this.
0174  */
0175 CELER_FUNCTION auto
0176 SBEnergyDistHelper::calc_max_xs(SBTables const& xs_params,
0177                                 ElementId element) const -> Xs
0178 {
0179     CELER_EXPECT(element);
0180     SBElementTableData const& el = xs_params.elements[element];
0181
0182     size_type const x_idx = calc_xs_.x_index();
0183     real_type const x_frac = calc_xs_.x_fraction();
0184     real_type result;
0185
0186     // Calc max xs
0187     auto get_value = [&xs_params, &el](size_type ix) -> real_type {
0188         // Index of the largest xs for exiting energy for the given
0189         // incident grid point
0190         size_type iy = xs_params.sizes[el.argmax[ix]];
0191         // Value of the maximum cross section
0192         return xs_params.reals[el.grid.at(ix, iy)];
0193     };
0194     result = (1 - x_frac) * get_value(x_idx) + x_frac * get_value(x_idx + 1);
0195
0196     CELER_ENSURE(result > 0);
0197     return Xs{result};
0198 }
0199
0200 //---------------------------------------------------------------------------//
0201 /*!
0202  * Construct a sampler for scaled exiting energy.
0203  */
0204 CELER_FUNCTION auto SBEnergyDistHelper::make_esq_sampler(
0205     real_type inc_energy, real_type min_gamma_energy) const -> ReciprocalSampler
0206 {
0207     CELER_EXPECT(min_gamma_energy > 0);
0208     return ReciprocalSampler(ipow<2>(min_gamma_energy) + dens_corr_,
0209                              ipow<2>(inc_energy) + dens_corr_);
0210 }
0211
0212 //---------------------------------------------------------------------------//
0213 }  // namespace celeritas