em/distribution/SBEnergyDistHelper.hh

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