model/detail/NuclearZoneBuilder.hh

0001 //----------------------------------*-C++-*----------------------------------//
0002 // Copyright 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/neutron/model/detail/NuclearZoneBuilder.hh
0007 //---------------------------------------------------------------------------//
0008 #pragma once
0009
0010 #include <cmath>
0011 #include <numeric>
0012 #include <vector>
0013
0014 #include "corecel/cont/Range.hh"
0015 #include "corecel/cont/Span.hh"
0016 #include "corecel/math/Integrator.hh"
0017 #include "celeritas/Types.hh"
0018 #include "celeritas/mat/IsotopeView.hh"
0019 #include "celeritas/neutron/data/NeutronInelasticData.hh"
0020 #include "celeritas/phys/AtomicNumber.hh"
0021
0022 #include "../CascadeOptions.hh"
0023
0024 namespace celeritas
0025 {
0026 namespace detail
0027 {
0028 //---------------------------------------------------------------------------//
0029 /*!
0030  * Construct NuclearZoneData for NeutronInelasticModel.
0031  */
0032 class NuclearZoneBuilder
0033 {
0034   public:
0035     //!@{
0036     //! \name Type aliases
0037     using AtomicMassNumber = AtomicNumber;
0038     using MevMass = units::MevMass;
0039     using Energy = units::MevEnergy;
0040     using ComponentVec = std::vector<ZoneComponent>;
0041     using Data = HostVal<NuclearZoneData>;
0042     //!@}
0043
0044   public:
0045     // Construct with cascade options and shared data
0046     inline NuclearZoneBuilder(CascadeOptions const& options,
0047                               NeutronInelasticScalars const& scalars,
0048                               Data* data);
0049
0050     // Construct nuclear zone data for a target (isotope)
0051     inline void operator()(IsotopeView const& target);
0052
0053   private:
0054     //// TYPES ////
0055
0056     struct ZoneDensity
0057     {
0058         real_type radius;
0059         real_type integral;
0060     };
0061     using VecZoneDensity = std::vector<ZoneDensity>;
0062
0063     //// DATA ////
0064
0065     // Cascade model configurations and nuclear structure parameters
0066     CascadeOptions const& options_;
0067
0068     real_type skin_depth_;
0069     MevMass neutron_mass_;
0070     MevMass proton_mass_;
0071
0072     CollectionBuilder<ZoneComponent> components_;
0073     CollectionBuilder<NuclearZones, MemSpace::host, IsotopeId> zones_;
0074
0075     //// HELPER FUNCTIONS ////
0076
0077     // Calculate components of nuclear zone properties
0078     inline ComponentVec calc_zone_components(IsotopeView const& target) const;
0079
0080     // Calculate zone densities for lightweight nuclei (A < 5)
0081     VecZoneDensity calc_zones_light(AtomicMassNumber a) const;
0082
0083     // Calculate zone densities for small nuclei (5 <= A < 12)
0084     VecZoneDensity calc_zones_small(AtomicMassNumber a) const;
0085
0086     // Calculate zone densities for heavier nuclei (A >= 12)
0087     VecZoneDensity calc_zones_heavy(AtomicMassNumber a) const;
0088
0089     // Calculate the nuclear radius
0090     inline real_type calc_nuclear_radius(AtomicMassNumber a) const;
0091 };
0092
0093 //---------------------------------------------------------------------------//
0094 // INLINE DEFINITIONS
0095 //---------------------------------------------------------------------------//
0096 /*!
0097  * Construct with cascade options and data.
0098  */
0099 NuclearZoneBuilder::NuclearZoneBuilder(CascadeOptions const& options,
0100                                        NeutronInelasticScalars const& scalars,
0101                                        Data* data)
0102     : options_(options)
0103     , skin_depth_(0.611207 * options.radius_scale)
0104     , neutron_mass_(scalars.neutron_mass)
0105     , proton_mass_(scalars.proton_mass)
0106     , components_(&data->components)
0107     , zones_(&data->zones)
0108 {
0109     CELER_EXPECT(options_);
0110     CELER_EXPECT(data);
0111 }
0112
0113 //---------------------------------------------------------------------------//
0114 /*!
0115  * Construct nuclear zone data for a single isotope.
0116  */
0117 void NuclearZoneBuilder::operator()(IsotopeView const& target)
0118 {
0119     auto comp = this->calc_zone_components(target);
0120
0121     NuclearZones nucl_zones;
0122     nucl_zones.zones = components_.insert_back(comp.begin(), comp.end());
0123     zones_.push_back(nucl_zones);
0124 }
0125
0126 //---------------------------------------------------------------------------//
0127 /*!
0128  * Calculate components of nuclear zone data.
0129  *
0130  * The nuclear zone radius, volume,
0131  * density, Fermi momentum and potential function as in G4NucleiModel and as
0132  * documented in section 24.2.3 of the Geant4 Physics Reference (release 11.2).
0133  *
0134  */
0135 auto NuclearZoneBuilder::calc_zone_components(IsotopeView const& target) const
0136     -> ComponentVec
0137 {
0138     using A = AtomicMassNumber;
0139
0140     // Calculate nuclear radius
0141     A const a = target.atomic_mass_number();
0142
0143     // Calculate per-zone radius and potential integral
0144     VecZoneDensity zone_dens;
0145     if (a < A{5})
0146     {
0147         zone_dens = this->calc_zones_light(a);
0148     }
0149     else if (a < A{12})
0150     {
0151         zone_dens = this->calc_zones_small(a);
0152     }
0153     else
0154     {
0155         zone_dens = this->calc_zones_heavy(a);
0156     }
0157     CELER_ASSERT(!zone_dens.empty());
0158
0159     // Fill the differential nuclear volume by each zone
0160     constexpr real_type four_thirds_pi = 4 * constants::pi / real_type{3};
0161
0162     ComponentVec components(zone_dens.size());
0163     real_type prev_volume{0};
0164     for (auto i : range(components.size()))
0165     {
0166         real_type volume = four_thirds_pi * ipow<3>(zone_dens[i].radius);
0167
0168         // Save radius and differential volume
0169         components[i].radius = zone_dens[i].radius;
0170         components[i].volume = volume - prev_volume;
0171         prev_volume = volume;
0172     }
0173
0174     // Fill the nuclear zone density, fermi momentum, and potential
0175     int num_protons = target.atomic_number().get();
0176     int num_nucleons[] = {num_protons, a.get() - num_protons};
0177     real_type mass[] = {proton_mass_.value(), neutron_mass_.value()};
0178     real_type dm[] = {target.proton_loss_energy().value(),
0179                       target.neutron_loss_energy().value()};
0180     static_assert(std::size(mass) == ZoneComponent::NucleonArray::size());
0181     static_assert(std::size(dm) == std::size(mass));
0182
0183     real_type const total_integral
0184         = std::accumulate(zone_dens.begin(),
0185                           zone_dens.end(),
0186                           real_type{0},
0187                           [](real_type sum, ZoneDensity const& zone) {
0188                               return sum + zone.integral;
0189                           });
0190
0191     for (auto i : range(components.size()))
0192     {
0193         for (auto ptype : range(ZoneComponent::NucleonArray::size()))
0194         {
0195             components[i].density[ptype]
0196                 = static_cast<real_type>(num_nucleons[ptype])
0197                   * zone_dens[i].integral
0198                   / (total_integral * components[i].volume);
0199             components[i].fermi_mom[ptype]
0200                 = options_.fermi_scale
0201                   * std::cbrt(components[i].density[ptype]);
0202             components[i].potential[ptype]
0203                 = ipow<2>(components[i].fermi_mom[ptype]) / (2 * mass[ptype])
0204                   + dm[ptype];
0205         }
0206     }
0207     return components;
0208 }
0209
0210 //---------------------------------------------------------------------------//
0211 /*!
0212  * Lightweight nuclei are treated as simple balls.
0213  */
0214 auto NuclearZoneBuilder::calc_zones_light(AtomicMassNumber a) const
0215     -> VecZoneDensity
0216 {
0217     CELER_EXPECT(a <= AtomicMassNumber{4});
0218     ZoneDensity result;
0219     result.radius = options_.radius_small
0220                     * (a == AtomicMassNumber{4} ? options_.radius_alpha : 1);
0221     result.integral = 1;
0222     return {result};
0223 }
0224
0225 //---------------------------------------------------------------------------//
0226 /*!
0227  * Small nuclei have a three-zone gaussian potential.
0228  */
0229 auto NuclearZoneBuilder::calc_zones_small(AtomicMassNumber a) const
0230     -> VecZoneDensity
0231 {
0232     real_type const nuclear_radius = this->calc_nuclear_radius(a);
0233     real_type gauss_radius = std::sqrt(
0234         ipow<2>(nuclear_radius) * (1 - 1 / static_cast<real_type>(a.get()))
0235         + real_type{6.4});
0236
0237     Integrator integrate_gauss{
0238         [](real_type r) { return ipow<2>(r) * std::exp(-ipow<2>(r)); }};
0239
0240     // Precompute y = sqrt(-log(alpha)) where alpha[3] = {0.7, 0.3, 0.01}
0241     real_type const y[] = {0.597223, 1.09726, 2.14597};
0242     VecZoneDensity result(std::size(y));
0243
0244     real_type ymin = 0;
0245     for (auto i : range(result.size()))
0246     {
0247         result[i].radius = gauss_radius * y[i];
0248         result[i].integral = ipow<3>(gauss_radius)
0249                              * integrate_gauss(ymin, y[i]);
0250         ymin = y[i];
0251     }
0252
0253     CELER_ENSURE(result.size() == 3);
0254     return result;
0255 }
0256
0257 //---------------------------------------------------------------------------//
0258 /*!
0259  * Heavy nuclei have a three- or six-zone Woods-Saxon potential.
0260  *
0261  * The Woods-Saxon potential, \f$ V(r) \f$,
0262  *
0263  * \f[
0264      V(r) = frac{V_{o}}{1 + e^{\frac{r - R}{a}}}
0265    \f]
0266  * is integrated numerically over the volume from \f$ r_{min} \f$ to \f$
0267  * r_{rmax} \f$, where \f$ V_{o}, R, a\f$ are the potential well depth, nuclear
0268  * radius, and surface thickness (skin depth), respectively.
0269  */
0270 auto NuclearZoneBuilder::calc_zones_heavy(AtomicMassNumber a) const
0271     -> VecZoneDensity
0272 {
0273     real_type const nuclear_radius = this->calc_nuclear_radius(a);
0274     real_type const skin_ratio = nuclear_radius / skin_depth_;
0275     real_type const skin_decay = std::exp(-skin_ratio);
0276     Integrator integrate_ws{[ws_shift = 2 * skin_ratio](real_type r) {
0277         return r * (r + ws_shift) / (1 + std::exp(r));
0278     }};
0279
0280     Span<real_type const> alpha;
0281     if (a < AtomicMassNumber{100})
0282     {
0283         static real_type const alpha_i[] = {0.7, 0.3, 0.01};
0284         alpha = make_span(alpha_i);
0285     }
0286     else
0287     {
0288         static real_type const alpha_h[] = {0.9, 0.6, 0.4, 0.2, 0.1, 0.05};
0289         alpha = make_span(alpha_h);
0290     }
0291
0292     VecZoneDensity result(alpha.size());
0293     real_type ymin = -skin_ratio;
0294     for (auto i : range(result.size()))
0295     {
0296         real_type y = std::log((1 + skin_decay) / alpha[i] - 1);
0297         result[i].radius = nuclear_radius + skin_depth_ * y;
0298
0299         result[i].integral
0300             = ipow<3>(skin_depth_)
0301               * (integrate_ws(ymin, y)
0302                  + ipow<2>(skin_ratio)
0303                        * std::log((1 + std::exp(-ymin)) / (1 + std::exp(-y))));
0304         ymin = y;
0305     }
0306
0307     return result;
0308 }
0309
0310 //---------------------------------------------------------------------------//
0311 /*!
0312  * Calculate the nuclear radius (R) computed from the atomic mass number (A).
0313  *
0314  * For \f$ A > 4 \f$, the nuclear radius with two parameters takes the form,
0315  * \f[
0316      R = [ 1.16 * A^{1/3} - 1.3456 / A^{1/3} ] \cdot R_{scale}
0317    \f]
0318  * where \f$ R_{scale} \f$ is a configurable parameter in [femtometer], while
0319  * \f$ R = 1.2 A^{1/3} \cdot R_{scale} \f$ (default) with a single parameter.
0320  */
0321 real_type NuclearZoneBuilder::calc_nuclear_radius(AtomicMassNumber a) const
0322 {
0323     CELER_EXPECT(a > AtomicMassNumber{4});
0324
0325     // Nuclear radius computed from A
0326     real_type cbrt_a = std::cbrt(static_cast<real_type>(a.get()));
0327     real_type par_a = (options_.use_two_params ? 1.16 : 1.2);
0328     real_type par_b = (options_.use_two_params ? -1.3456 : 0);
0329     return options_.radius_scale * (par_a * cbrt_a + par_b / cbrt_a);
0330 }
0331
0332 //---------------------------------------------------------------------------//
0333 }  // namespace detail
0334 }  // namespace celeritas