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