diff --git a/doc/_static/zotero.bib b/doc/_static/zotero.bib index 36c84fa9a9..bc5af07237 100644 --- a/doc/_static/zotero.bib +++ b/doc/_static/zotero.bib @@ -3747,32 +3747,64 @@ @article{bom-experimentaldt-2005 } @article{yamashita-muonicspin-2022, - title = {Roles of resonant muonic molecule in new kinetics model and muon catalyzed fusion in compressed gas}, - author = {Yamashita, Takuma and Kino, Yasushi and Okutsu, Kenichi and Okada, Shinji and Sato, Motoyasu}, - year = 2022, - month = apr, - volume = {12}, - number = {1}, - journal = {Scientific Reports}, - shortjournal = {Sci Rep}, - issn = {2045-2322}, - doi = {10.1038/s41598-022-09487-0}, - url = {https://www.nature.com/articles/s41598-022-09487-0}, - pages = {6393} + title = {Roles of resonant muonic molecule in new kinetics model and muon catalyzed fusion in compressed gas}, + author = {Yamashita, Takuma and Kino, Yasushi and Okutsu, Kenichi and Okada, Shinji and Sato, Motoyasu}, + year = 2022, + month = apr, + volume = {12}, + number = {1}, + journal = {Scientific Reports}, + shortjournal = {Sci Rep}, + issn = {2045-2322}, + doi = {10.1038/s41598-022-09487-0}, + url = {https://www.nature.com/articles/s41598-022-09487-0}, + pages = {6393} } @article{faifman-mucfformation-1996, - title = {Quadrupole corrections to matrix elements of transitions in resonant reactions of muonic molecule formation}, - author = {Faifman, M. P. and Strizh, T. A. and Armour, E. A. G. and Harston, M. R.}, - year = 1996, - month = dec, - volume = {101-102}, - rights = {http://www.springer.com/tdm}, - issn = {0304-3834, 1572-9540}, - url = {http://link.springer.com/10.1007/BF02227621}, - doi = {10.1007/BF02227621}, - pages = {179--189}, - number = {1}, - journaltitle = {Hyperfine Interactions}, - shortjournal = {Hyperfine Interact}, + title = {Quadrupole corrections to matrix elements of transitions in resonant reactions of muonic molecule formation}, + author = {Faifman, M. P. and Strizh, T. A. and Armour, E. A. G. and Harston, M. R.}, + year = 1996, + month = dec, + volume = {101-102}, + rights = {http://www.springer.com/tdm}, + issn = {0304-3834, 1572-9540}, + url = {http://link.springer.com/10.1007/BF02227621}, + doi = {10.1007/BF02227621}, + pages = {179--189}, + number = {1}, + journaltitle = {Hyperfine Interactions}, + shortjournal = {Hyperfine Interact} +} + +@article{balin-mucf-2011, + title = {High Precision Study of Muon Catalyzed Fusion in {{D2}} and {{HD}} Gas}, + author = {Balin, D. V. and Ganzha, V. A. and Kozlov, S. M. and Maev, E. M. and Petrov, G. E. and Soroka, M. A. and Schapkin, G. N. and Semenchuk, G. G. and Trofimov, V. A. and Vasiliev, A. A. and Vorobyov, A. A. and Voropaev, N. I. and Petitjean, C. and Gartner, B. and Lauss, B. and Marton, J. and Zmeskal, J. and Case, T. and Crowe, K. M. and Kammel, P. and Hartmann, F. J. and Faifman, M. P.}, + year = 2011, + month = mar, + journal = {Physics of Particles and Nuclei}, + volume = {42}, + number = {2}, + pages = {185--214}, + issn = {1063-7796, 1531-8559}, + doi = {10.1134/S106377961102002X}, + urldate = {2026-02-20}, + copyright = {http://www.springer.com/tdm}, + langid = {english} +} + +@article{bogdanova-mucf-2009, + title = {{Experimental Investigation of Muon-Catalyzed t + t Fusion}}, + author = {Bogdanova, L. N. and Bom, V. R. and Demin, A. M. and Demin, D. L. and Van Eijk, C. W. E. and Filchagin, S. V. and Filchenkov, V. V. and Grafov, N. N. and Grishechkin, S. K. and Gritsaj, K. I. and Konin, A. D. and Kuryakin, A. V. and Medved', S. V. and Musyaev, R. K. and Rudenko, A. I. and Tumkin, D. P. and Vinogradov, {\relax Yu}. I. and Yukhimchuk, A. A. and Yukhimchuk, S. A. and Zinov, V. G. and Zlatoustovskii, S. V.}, + year = 2009, + month = feb, + journal = {Journal of Experimental and Theoretical Physics}, + volume = {108}, + number = {2}, + pages = {216--225}, + issn = {1063-7761, 1090-6509}, + doi = {10.1134/S1063776109020034}, + urldate = {2026-02-20}, + copyright = {http://www.springer.com/tdm}, + langid = {english} } diff --git a/doc/implementation/mucf-physics.rst b/doc/implementation/mucf-physics.rst index 819c263f85..390d3fcd40 100644 --- a/doc/implementation/mucf-physics.rst +++ b/doc/implementation/mucf-physics.rst @@ -133,14 +133,18 @@ host/device data via :cpp:class:`celeritas::detail::MucfMaterialInserter`. .. doxygenclass:: celeritas::detail::MucfMaterialInserter The main cycle is managed by the model's -:cpp:class:`celeritas::DTMixMucfExecutor`. The muonic atom selection is handled -by detail classes, while the Interactors are reserved for sampling final states -of the outgoing secondaries. +:cpp:class:`celeritas::DTMixMucfExecutor`. The muonic atom and molecule +formation are handled by detail classes, while the Interactors are reserved for +sampling final states of the outgoing secondaries. .. note:: Only reactive channels are implemented. .. doxygenclass:: celeritas::detail::MuonicAtomSelector .. doxygenclass:: celeritas::detail::MuonicAtomSpinSelector +.. doxygenclass:: celeritas::detail::MuonicMoleculeSelector +.. doxygenclass:: celeritas::detail::DDChannelSelector +.. doxygenclass:: celeritas::detail::DTChannelSelector +.. doxygenclass:: celeritas::detail::TTChannelSelector .. doxygenclass:: celeritas::DDMucfInteractor .. doxygenclass:: celeritas::DTMucfInteractor .. doxygenclass:: celeritas::TTMucfInteractor diff --git a/src/celeritas/inp/MucfPhysics.cc b/src/celeritas/inp/MucfPhysics.cc index 22cf2a963e..d5ede1fd98 100644 --- a/src/celeritas/inp/MucfPhysics.cc +++ b/src/celeritas/inp/MucfPhysics.cc @@ -503,7 +503,8 @@ MucfCycleRate hd_1_cycle_data() /*! * Muon-catalyzed fusion cycle rate data for \em tt fusion with F = 1/2 state. * - * The tt fusion cycle rate data is currently a constant value at 2.8e6 1/s. + * The tt fusion cycle rate data is currently a constant value at 2.8e6 1/s + * citet{bogdanova-mucf-2009, https://doi.org/10.1134/S1063776109020034} . */ MucfCycleRate tt_1_over_2_cycle_data() { diff --git a/src/celeritas/mucf/Types.hh b/src/celeritas/mucf/Types.hh index 23a02637c0..f23f7e3b80 100644 --- a/src/celeritas/mucf/Types.hh +++ b/src/celeritas/mucf/Types.hh @@ -7,6 +7,7 @@ #pragma once #include "corecel/OpaqueId.hh" +#include "celeritas/Quantities.hh" namespace celeritas { @@ -55,7 +56,23 @@ enum class MucfMuonicMolecule //---------------------------------------------------------------------------// //! Opaque index of a muCF material component -using MuCfMatId = OpaqueId; +using MucfMatId = OpaqueId; + +//---------------------------------------------------------------------------// +// REACTIVE SPIN STATES +//---------------------------------------------------------------------------// + +//! State F = 1/2 (dd and tt) +inline constexpr units::HalfSpinInt spin_one_half{1}; + +//! State F = 3/2 (dd) +inline constexpr units::HalfSpinInt spin_three_halves{3}; + +//! State F = 0 (dt) +inline constexpr units::HalfSpinInt spin_zero{0}; + +//! State F = 1 (dt) +inline constexpr units::HalfSpinInt spin_one{2}; //---------------------------------------------------------------------------// } // namespace celeritas diff --git a/src/celeritas/mucf/data/DTMixMucfData.hh b/src/celeritas/mucf/data/DTMixMucfData.hh index ff6aeb7c66..255effccf6 100644 --- a/src/celeritas/mucf/data/DTMixMucfData.hh +++ b/src/celeritas/mucf/data/DTMixMucfData.hh @@ -106,9 +106,9 @@ struct DTMixMucfData template using Items = Collection; template - using MaterialItems = Collection; + using MaterialItems = Collection; using GridRecord = NonuniformGridRecord; - using CycleTimesArray = EnumArray>; + using CycleRatesArray = EnumArray>; using MaterialFractionsArray = EnumArray; //! Particles @@ -122,12 +122,12 @@ struct DTMixMucfData //!@{ //! Material-dependent data calculated at model construction - //! \c PhysMatId indexed by \c MuCfMatId + //! \c PhysMatId indexed by \c MucfMatId MaterialItems mucfmatid_to_matid; //! Isotopic fractions per material: [mat_comp_id][isotope] MaterialItems isotopic_fractions; - //! Cycle times per material: [mat_comp_id][muonic_molecule][spin_index] - MaterialItems cycle_times; //!< In [s] + //! Cycle rates per material: [mat_comp_id][muonic_molecule][spin_index] + MaterialItems cycle_rates; //!< In [1/s] //! \todo Add mean atom spin flip times //! \todo Add mean atom transfer times //!@} @@ -137,7 +137,7 @@ struct DTMixMucfData { return particle_ids && particle_masses && muon_energy_cdf && !reals.empty() && !mucfmatid_to_matid.empty() - && !isotopic_fractions.empty() && !cycle_times.empty(); + && !isotopic_fractions.empty() && !cycle_rates.empty(); } //! Assign from another set of data @@ -153,7 +153,7 @@ struct DTMixMucfData this->reals = other.reals; this->mucfmatid_to_matid = other.mucfmatid_to_matid; this->isotopic_fractions = other.isotopic_fractions; - this->cycle_times = other.cycle_times; + this->cycle_rates = other.cycle_rates; return *this; } diff --git a/src/celeritas/mucf/executor/DTMixMucfExecutor.hh b/src/celeritas/mucf/executor/DTMixMucfExecutor.hh index d2326bf93d..000f7ed235 100644 --- a/src/celeritas/mucf/executor/DTMixMucfExecutor.hh +++ b/src/celeritas/mucf/executor/DTMixMucfExecutor.hh @@ -16,10 +16,9 @@ #include "detail/DDChannelSelector.hh" #include "detail/DTChannelSelector.hh" -#include "detail/DTMixMuonicMoleculeSelector.hh" #include "detail/MuonicAtomSelector.hh" #include "detail/MuonicAtomSpinSelector.hh" -#include "detail/MuonicMoleculeSpinSelector.hh" +#include "detail/MuonicMoleculeSelector.hh" #include "detail/TTChannelSelector.hh" namespace celeritas @@ -44,39 +43,36 @@ DTMixMucfExecutor::operator()(celeritas::CoreTrackView const& track) auto elcomp_id = phys_step_view.element(); CELER_ASSERT(elcomp_id); - auto element = track.material().material_record().element_record(elcomp_id); + auto const& mat_record = track.material().material_record(); + auto element = mat_record.element_record(elcomp_id); CELER_ASSERT(element.atomic_number() == AtomicNumber{1}); // Must be H - auto rng = track.rng(); - - // Muon decay may compete against other "actions" in this executor - real_type const decay_len{}; //! \todo Set muon decay interaction length + //! \todo Make sure that at this point we selected d or t already // Find muCF material ID from PhysMatId // Make this a View if ever used beyond this executor - auto find = [&](PhysMatId matid) -> MuCfMatId { + auto find = [&](PhysMatId matid) -> MucfMatId { CELER_EXPECT(matid); for (auto i : range(data.mucfmatid_to_matid.size())) { - if (auto const comp_id = MuCfMatId{i}; + if (auto const comp_id = MucfMatId{i}; data.mucfmatid_to_matid[comp_id] == matid) { return comp_id; } } // MuCF material ID not found - return MuCfMatId{}; + return MucfMatId{}; }; auto const mucf_matid = find(track.material().material_id()); CELER_ASSERT(mucf_matid); - // Form d or t muonic atom - detail::MuonicAtomSelector form_atom( - data.isotopic_fractions[mucf_matid][MucfIsotope::deuterium]); - auto muonic_atom = form_atom(rng); + auto rng = track.rng(); - // Select atom spin via a helper class - detail::MuonicAtomSpinSelector select_atom_spin(muonic_atom); + // Form d or t muonic atom + auto muonic_atom = detail::MuonicAtomSelector( + data.isotopic_fractions[mucf_matid][MucfIsotope::deuterium])(rng); + auto atom_spin = detail::MuonicAtomSpinSelector(muonic_atom)(rng); // { // Competing at-rest processes which add to the total track time @@ -85,31 +81,13 @@ DTMixMucfExecutor::operator()(celeritas::CoreTrackView const& track) // } // Form dd, dt, or tt muonic molecule - detail::DTMixMuonicMoleculeSelector form_muonic_molecule; - auto muonic_molecule = form_muonic_molecule(rng); - - // Select molecule spin - detail::MuonicMoleculeSpinSelector select_molecule_spin(muonic_molecule); - auto const molecule_spin = select_molecule_spin(rng); - - // Load cycle time for the selected molecule - auto const cycle_time - = data.cycle_times[mucf_matid][muonic_molecule][molecule_spin]; - CELER_ASSERT(cycle_time > 0); - - // Check if muon decays before fusion happens - real_type const mucf_len = cycle_time * track.sim().step_length(); - if (decay_len < mucf_len) - { - // Muon decays and halts the interaction - //! \todo Update track time and return muon decay interactor - } + auto [muonic_molecule, cycle_time] = detail::MuonicMoleculeSelector( + muonic_atom, atom_spin, data.cycle_rates[mucf_matid])(rng); - //! \todo Correct track time update? Or should be done in Interactors? + // Update track time according to the sampled cycle time track.sim().add_time(cycle_time); // Fuse molecule and generate secondaries - //! \todo Maybe move the channel selectors into the interactors auto allocate_secondaries = phys_step_view.make_secondary_allocator(); Interaction result; switch (muonic_molecule) @@ -117,7 +95,9 @@ DTMixMucfExecutor::operator()(celeritas::CoreTrackView const& track) case MucfMuonicMolecule::deuterium_deuterium: { // Return DD interaction DDMucfInteractor interact( - data, detail::DDChannelSelector()(rng), allocate_secondaries); + data, + detail::DDChannelSelector(mat_record.temperature())(rng), + allocate_secondaries); result = interact(rng); break; } diff --git a/src/celeritas/mucf/executor/detail/DDChannelSelector.hh b/src/celeritas/mucf/executor/detail/DDChannelSelector.hh index 05394c8019..11b896c23c 100644 --- a/src/celeritas/mucf/executor/detail/DDChannelSelector.hh +++ b/src/celeritas/mucf/executor/detail/DDChannelSelector.hh @@ -6,6 +6,7 @@ //---------------------------------------------------------------------------// #pragma once +#include "corecel/random/distribution/BernoulliDistribution.hh" #include "celeritas/mucf/interactor/DDMucfInteractor.hh" namespace celeritas @@ -16,8 +17,13 @@ namespace detail /*! * Select final channel for muonic dd molecules. * - * This selection already accounts for sticking, as that is one of the possible - * outcomes. + * The material temperature determines the probability + * of the outcome of the fusion ending in the \f$ ^3\text{He} \f$ channels + * versus the tritium channel. If the outcome is a \f$ ^3\text{He} \f$ channel, + * a constant sticking fraction of 12.2% is used to define if sticking occurs + * \citet{balin-mucf-2011, https://doi.org/10.1134/S106377961102002X} . + * + * \todo Update I/O with user-defined sticking fractions and branching ratios. */ class DDChannelSelector { @@ -27,41 +33,75 @@ class DDChannelSelector using Channel = DDMucfInteractor::Channel; //!@} - public: - //! Construct with args; \todo Update documentation - inline CELER_FUNCTION DDChannelSelector(/* args */); + //! Construct with material temperature + inline CELER_FUNCTION DDChannelSelector(real_type material_temperature); // Select fusion channel to be used by the interactor template - inline CELER_FUNCTION Channel operator()(Engine&); + inline CELER_FUNCTION Channel operator()(Engine& rng); + + private: + real_type he3_probability_; + + // Constant sticking fraction between the two 3He channels + inline CELER_FUNCTION static constexpr real_type sticking_fraction() + { + return 0.122; + } }; //---------------------------------------------------------------------------// // INLINE DEFINITIONS //---------------------------------------------------------------------------// /*! - * Construct with args. + * Construct with material temperature. * - * \todo Update documentation + * The temperature is used to calculate the branching ratio between either of + * the \f$ ^3\text{He} \f$ channels and the tritium channel. */ -CELER_FUNCTION DDChannelSelector::DDChannelSelector(/* args */) +CELER_FUNCTION +DDChannelSelector::DDChannelSelector(real_type material_temperature) { - CELER_NOT_IMPLEMENTED("Mucf dd fusion channel selection"); + CELER_EXPECT(material_temperature > 0); + + real_type branching_ratio{0}; + if (material_temperature <= 50) + { + branching_ratio = 1; + } + else if (material_temperature <= 100) + { + branching_ratio = 1.0088 * (material_temperature - 50); + } + else + { + branching_ratio = 1.44; + } + CELER_ASSERT(branching_ratio > 0); + + he3_probability_ = branching_ratio / (branching_ratio + 1); + CELER_ENSURE(he3_probability_ > 0 && he3_probability_ <= 1); } //---------------------------------------------------------------------------// /*! - * Return a sampled channel to be used as input in the dd muCF interactor. + * Return a selected fusion channel for the \f$ (dd)_\mu \f$ muonic molecule. * * \sa celeritas::DDMucfInteractor */ template -CELER_FUNCTION DDChannelSelector::Channel DDChannelSelector::operator()(Engine&) +CELER_FUNCTION DDChannelSelector::Channel +DDChannelSelector::operator()(Engine& rng) { Channel result{Channel::size_}; - //! \todo Implement - // Final channel selection already takes into account sticking. + result = BernoulliDistribution(he3_probability_)(rng) + // Select between the two 3He channels + ? (BernoulliDistribution(this->sticking_fraction())(rng) + ? Channel::muonichelium3_neutron + : Channel::helium3_muon_neutron) + // Select tritium channel + : Channel::tritium_muon_proton; CELER_ENSURE(result < Channel::size_); return result; diff --git a/src/celeritas/mucf/executor/detail/DTChannelSelector.hh b/src/celeritas/mucf/executor/detail/DTChannelSelector.hh index 47e097cdfd..b890d931f7 100644 --- a/src/celeritas/mucf/executor/detail/DTChannelSelector.hh +++ b/src/celeritas/mucf/executor/detail/DTChannelSelector.hh @@ -6,6 +6,7 @@ //---------------------------------------------------------------------------// #pragma once +#include "corecel/random/distribution/BernoulliDistribution.hh" #include "celeritas/mucf/interactor/DTMucfInteractor.hh" namespace celeritas @@ -16,8 +17,11 @@ namespace detail /*! * Select final channel for muonic dt molecules. * - * This selection already accounts for sticking, as that is one of the possible - * outcomes. + * The selection is based on a constant sticking fraction from + * \citet{kamimura-mucf-2023, https://doi.org/10.1103/PhysRevC.107.034607} + * in which ~0.8% of the time the muonic alpha channel is selected. + * + * \todo Update I/O with user-defined sticking fractions. */ class DTChannelSelector { @@ -27,41 +31,38 @@ class DTChannelSelector using Channel = DTMucfInteractor::Channel; //!@} - public: - //! Construct with args; \todo Update documentation - inline CELER_FUNCTION DTChannelSelector(/* args */); + //! Default constructor + inline CELER_FUNCTION DTChannelSelector() = default; // Select fusion channel to be used by the interactor template - inline CELER_FUNCTION Channel operator()(Engine&); + inline CELER_FUNCTION Channel operator()(Engine& rng); + + private: + // Constant sticking fraction of dt fusion + inline CELER_FUNCTION static constexpr real_type sticking_fraction() + { + return 0.00857; + } }; //---------------------------------------------------------------------------// // INLINE DEFINITIONS //---------------------------------------------------------------------------// /*! - * Construct with args. - * - * \todo Update documentation - */ -CELER_FUNCTION DTChannelSelector::DTChannelSelector(/* args */) -{ - CELER_NOT_IMPLEMENTED("Mucf dt fusion channel selection"); -} - -//---------------------------------------------------------------------------// -/*! - * Return a sampled channel to be used as input in the dt muCF interactor. + * Return a selected fusion channel for the \f$ (dt)_\mu \f$ muonic molecule. * * \sa celeritas::DTMucfInteractor */ template -CELER_FUNCTION DTChannelSelector::Channel DTChannelSelector::operator()(Engine&) +CELER_FUNCTION DTChannelSelector::Channel +DTChannelSelector::operator()(Engine& rng) { Channel result{Channel::size_}; - //! \todo Implement - // Final channel selection already takes into account sticking. + result = (BernoulliDistribution(this->sticking_fraction())(rng)) + ? Channel::muonicalpha_neutron + : Channel::alpha_muon_neutron; CELER_ENSURE(result < Channel::size_); return result; diff --git a/src/celeritas/mucf/executor/detail/DTMixMuonicMoleculeSelector.hh b/src/celeritas/mucf/executor/detail/DTMixMuonicMoleculeSelector.hh deleted file mode 100644 index cc4c7bd30b..0000000000 --- a/src/celeritas/mucf/executor/detail/DTMixMuonicMoleculeSelector.hh +++ /dev/null @@ -1,67 +0,0 @@ -//------------------------------- -*- C++ -*- -------------------------------// -// Copyright Celeritas contributors: see top-level COPYRIGHT file for details -// SPDX-License-Identifier: (Apache-2.0 OR MIT) -//---------------------------------------------------------------------------// -//! \file celeritas/mucf/executor/detail/DTMixMuonicMoleculeSelector.hh -//---------------------------------------------------------------------------// -#pragma once - -#include "celeritas/Types.hh" - -namespace celeritas -{ -namespace detail -{ -//---------------------------------------------------------------------------// -/*! - * Select a muonic molecule by calculating the interaction lengths of the - * possible molecule formations. - * - * This is the equivalent of Geant4's - * \c G4VRestProcess::AtRestGetPhysicalInteractionLength - */ -class DTMixMuonicMoleculeSelector -{ - public: - //! Construct with args; \todo Update documentation - inline CELER_FUNCTION DTMixMuonicMoleculeSelector(/* args */); - - // Select muonic molecule - template - inline CELER_FUNCTION MucfMuonicMolecule operator()(Engine&); -}; - -//---------------------------------------------------------------------------// -// INLINE DEFINITIONS -//---------------------------------------------------------------------------// -/*! - * Construct with args. - * - * \todo Update documentation - */ -CELER_FUNCTION -DTMixMuonicMoleculeSelector::DTMixMuonicMoleculeSelector(/* args */) -{ - //! \todo Implement - CELER_NOT_IMPLEMENTED("Mucf muonic molecule selection"); -} - -//---------------------------------------------------------------------------// -/*! - * Return selected muonic molecule. - */ -template -CELER_FUNCTION MucfMuonicMolecule -DTMixMuonicMoleculeSelector::operator()(Engine&) -{ - MucfMuonicMolecule result{MucfMuonicMolecule::size_}; - - //! \todo Implement - - CELER_ENSURE(result < MucfMuonicMolecule::size_); - return result; -} - -//---------------------------------------------------------------------------// -} // namespace detail -} // namespace celeritas diff --git a/src/celeritas/mucf/executor/detail/MuonicMoleculeSelector.hh b/src/celeritas/mucf/executor/detail/MuonicMoleculeSelector.hh new file mode 100644 index 0000000000..c5e70045fc --- /dev/null +++ b/src/celeritas/mucf/executor/detail/MuonicMoleculeSelector.hh @@ -0,0 +1,189 @@ +//------------------------------- -*- C++ -*- -------------------------------// +// Copyright Celeritas contributors: see top-level COPYRIGHT file for details +// SPDX-License-Identifier: (Apache-2.0 OR MIT) +//---------------------------------------------------------------------------// +//! \file celeritas/mucf/executor/detail/MuonicMoleculeSelector.hh +//---------------------------------------------------------------------------// +#pragma once + +#include "corecel/cont/EnumArray.hh" +#include "corecel/cont/Range.hh" +#include "corecel/math/NumericLimits.hh" +#include "corecel/random/distribution/ExponentialDistribution.hh" +#include "celeritas/Quantities.hh" +#include "celeritas/Types.hh" +#include "celeritas/mucf/Types.hh" + +namespace celeritas +{ +namespace detail +{ +//---------------------------------------------------------------------------// +/*! + * Select a muonic molecule and its final cycle time. + * + * The muonic molecule is selected by sampling the shortest fusion cycle time + * \f[ + * T = -\ln(r) \times \tau_\text{cycle}, + * \f] + * + * where \f$ \tau_\text{cycle} \f$ is the cycle time a given molecule + spin + * calculated from material data. + */ +class MuonicMoleculeSelector +{ + public: + //!@{ + //! \name Type aliases + using CycleRatesArray = EnumArray>; + using HalfSpinInt = units::HalfSpinInt; + //!@} + + // Result of molecule and cycle time selection + struct Result + { + MucfMuonicMolecule molecule{MucfMuonicMolecule::size_}; + real_type cycle_time{numeric_limits::max()}; + + //! Check whether the data are assigned + explicit CELER_FUNCTION operator bool() const + { + return molecule < MucfMuonicMolecule::size_ && cycle_time > 0 + && cycle_time < numeric_limits::max(); + } + }; + + //! Construct with muonic atom and material information + inline CELER_FUNCTION MuonicMoleculeSelector(MucfMuonicAtom atom, + HalfSpinInt spin, + CycleRatesArray cycle_rates); + + // Select muonic molecule and sample its final cycle time + template + inline CELER_FUNCTION Result operator()(Engine& rng); + + private: + MucfMuonicAtom atom_; + size_type cycle_rate_index_; + CycleRatesArray cycle_rates_; + + // Sample the final cycle time for a given molecule + template + inline CELER_FUNCTION real_type sample_exp_time(MucfMuonicMolecule molecule, + Engine& rng); +}; + +//---------------------------------------------------------------------------// +// INLINE DEFINITIONS +//---------------------------------------------------------------------------// +/*! + * Construct with muonic atom and cycle rate information. + * + * The correct cycle rate array index is determined at construction by the + * combination of atom+spin. Given the limited number of states, the position + * in the index is set in \c MucfMaterialInserter::calc_[dd|dt|tt]_cycle + * manually. + * + * \note Current implementation has no safe method to access cycle rate data + * for each spin state. + */ +CELER_FUNCTION +MuonicMoleculeSelector::MuonicMoleculeSelector(MucfMuonicAtom atom, + HalfSpinInt spin, + CycleRatesArray cycle_rates) + : atom_(atom), cycle_rates_(cycle_rates) +{ + CELER_EXPECT(atom < MucfMuonicAtom::size_); + + // Check that the spin value is valid for the given atom type and set the + // cycle rate array index accordingly + switch (atom_) + { + case MucfMuonicAtom::deuterium: { + CELER_EXPECT(spin == spin_one_half || spin == spin_three_halves); + // DD fusion + // F = 1/2 and F = 3/2 correspond to indices 0 and 1, respectively + cycle_rate_index_ = (spin == spin_one_half) ? 0 : 1; + break; + } + case MucfMuonicAtom::tritium: { + CELER_EXPECT(spin == spin_zero || spin == spin_one + || spin == spin_one_half); + // DT: F = 0 and F = 1 correspond to indices 0 and 1, respectively + // TT: F = 1/2 corresponds to index 0 + cycle_rate_index_ + = (spin == spin_zero || spin == spin_one_half) ? 0 : 1; + break; + } + default: + CELER_ASSERT_UNREACHABLE(); + } +} + +//---------------------------------------------------------------------------// +/*! + * Return a muonic molecule and its sampled cycle time. + */ +template +CELER_FUNCTION MuonicMoleculeSelector::Result +MuonicMoleculeSelector::operator()(Engine& rng) +{ + using MMM = MucfMuonicMolecule; + using CycleTimeArray = EnumArray; + + CycleTimeArray sampled_times; + auto const inf = numeric_limits::max(); + if (atom_ == MucfMuonicAtom::deuterium) + { + // DD fusion is only triggered by a muonic deuterium + sampled_times[MMM::deuterium_deuterium] + = this->sample_exp_time(MMM::deuterium_deuterium, rng); + sampled_times[MMM::deuterium_tritium] = inf; + sampled_times[MMM::tritium_tritium] = inf; + } + else + { + // DT and TT fusions are triggered by a muonic tritium + sampled_times[MMM::deuterium_deuterium] = inf; + sampled_times[MMM::deuterium_tritium] + = this->sample_exp_time(MMM::deuterium_tritium, rng); + sampled_times[MMM::tritium_tritium] + = this->sample_exp_time(MMM::tritium_tritium, rng); + } + + Result result; + for (auto mol : range(MMM::size_)) + { + if (sampled_times[mol] < result.cycle_time) + { + result.cycle_time = sampled_times[mol]; + result.molecule = mol; + } + } + + CELER_ENSURE(result); + return result; +} + +//---------------------------------------------------------------------------// +/*! + * Sample the final exponential distribution for a given cycle rate using + * \f$ T = -\ln(r) \times \tau_\text{cycle} \f$ . If the cached rate is + * zero, the returned time is set to infinity, which effectively removes + * that molecule as a possible selection. + */ +template +CELER_FUNCTION real_type MuonicMoleculeSelector::sample_exp_time( + MucfMuonicMolecule molecule, Engine& rng) +{ + CELER_EXPECT(molecule < MucfMuonicMolecule::size_); + auto const rate = cycle_rates_[molecule][cycle_rate_index_]; + + // Return an infinite cycle time when the cached rate is zero + return (rate > 0) ? ExponentialDistribution(rate)(rng) + : numeric_limits::max(); +} + +//---------------------------------------------------------------------------// +} // namespace detail +} // namespace celeritas diff --git a/src/celeritas/mucf/executor/detail/MuonicMoleculeSpinSelector.hh b/src/celeritas/mucf/executor/detail/MuonicMoleculeSpinSelector.hh deleted file mode 100644 index 67ae8d4e29..0000000000 --- a/src/celeritas/mucf/executor/detail/MuonicMoleculeSpinSelector.hh +++ /dev/null @@ -1,64 +0,0 @@ -//------------------------------- -*- C++ -*- -------------------------------// -// Copyright Celeritas contributors: see top-level COPYRIGHT file for details -// SPDX-License-Identifier: (Apache-2.0 OR MIT) -//---------------------------------------------------------------------------// -//! \file celeritas/mucf/executor/detail/MuonicMoleculeSpinSelector.hh -//---------------------------------------------------------------------------// -#pragma once - -#include "corecel/Assert.hh" -#include "celeritas/Types.hh" - -namespace celeritas -{ -namespace detail -{ -//---------------------------------------------------------------------------// -/*! - * Select muonic molecule spin, in units of \f$ \frac{\hbar}{2} \f$. - * - * Applicable to \f$ (dd)_\mu \f$, * \f$ (dt)_\mu \f$ , and \f$ (tt)_\mu \f$. - */ -class MuonicMoleculeSpinSelector -{ - public: - // Construct with muonic molecule type - inline CELER_FUNCTION - MuonicMoleculeSpinSelector(MucfMuonicMolecule molecule); - - // Sample and return a spin value in units of hbar / 2 - template - inline CELER_FUNCTION size_type operator()(Engine&); - - private: - MucfMuonicMolecule molecule_; -}; - -//---------------------------------------------------------------------------// -// INLINE DEFINITIONS -//---------------------------------------------------------------------------// -/*! - * Construct with muonic molecule. - */ -CELER_FUNCTION MuonicMoleculeSpinSelector::MuonicMoleculeSpinSelector( - MucfMuonicMolecule molecule) - : molecule_(molecule) -{ - CELER_EXPECT(molecule_ < MucfMuonicMolecule::size_); - CELER_NOT_IMPLEMENTED("Mucf muonic molecule spin selection"); -} - -//---------------------------------------------------------------------------// -/*! - * Return selected spin, in units of \f$ \hbar / 2 \f$. - */ -template -CELER_FUNCTION size_type MuonicMoleculeSpinSelector::operator()(Engine&) -{ - //! \todo switch on molecule_ - CELER_ASSERT_UNREACHABLE(); -} - -//---------------------------------------------------------------------------// -} // namespace detail -} // namespace celeritas diff --git a/src/celeritas/mucf/executor/detail/TTChannelSelector.hh b/src/celeritas/mucf/executor/detail/TTChannelSelector.hh index dd8372479a..3716d073b4 100644 --- a/src/celeritas/mucf/executor/detail/TTChannelSelector.hh +++ b/src/celeritas/mucf/executor/detail/TTChannelSelector.hh @@ -6,6 +6,7 @@ //---------------------------------------------------------------------------// #pragma once +#include "corecel/random/distribution/BernoulliDistribution.hh" #include "celeritas/mucf/interactor/TTMucfInteractor.hh" namespace celeritas @@ -16,8 +17,11 @@ namespace detail /*! * Select final channel for muonic tt molecules. * - * This selection already accounts for sticking, as that is one of the possible - * outcomes. + * The selection is based on a constant sticking fraction from + * \citet{bogdanova-mucf-2009, https://doi.org/10.1134/S1063776109020034} , + * in which ~14% of the time the muonic alpha channel is selected. + * + * \todo Update I/O with user-defined sticking fractions. */ class TTChannelSelector { @@ -27,42 +31,38 @@ class TTChannelSelector using Channel = TTMucfInteractor::Channel; //!@} - public: - //! Construct with args; \todo Update documentation - inline CELER_FUNCTION TTChannelSelector(/* args */); + //! Default constructor + inline CELER_FUNCTION TTChannelSelector() = default; // Select fusion channel to be used by the interactor template - inline CELER_FUNCTION Channel operator()(Engine&); + inline CELER_FUNCTION Channel operator()(Engine& rng); + + private: + // Constant sticking fraction of tt fusion + inline CELER_FUNCTION real_type static constexpr sticking_fraction() + { + return 0.14; + } }; //---------------------------------------------------------------------------// // INLINE DEFINITIONS //---------------------------------------------------------------------------// /*! - * Construct with args. - * - * \todo Update documentation - */ -CELER_FUNCTION TTChannelSelector::TTChannelSelector(/* args */) -{ - //! \todo Implement - CELER_NOT_IMPLEMENTED("Mucf tt fusion channel selection"); -} - -//---------------------------------------------------------------------------// -/*! - * Return a sampled channel to be used as input in the tt muCF interactor. + * Return a selected fusion channel for the \f$ (tt)_\mu \f$ muonic molecule. * * \sa celeritas::TTMucfInteractor */ template -CELER_FUNCTION TTChannelSelector::Channel TTChannelSelector::operator()(Engine&) +CELER_FUNCTION TTChannelSelector::Channel +TTChannelSelector::operator()(Engine& rng) { Channel result{Channel::size_}; - //! \todo Implement - // Final channel selection already takes into account sticking. + result = (BernoulliDistribution(this->sticking_fraction())(rng)) + ? Channel::muonicalpha_neutron_neutron + : Channel::alpha_muon_neutron_neutron; CELER_ENSURE(result < Channel::size_); return result; diff --git a/src/celeritas/mucf/interactor/DTMucfInteractor.hh b/src/celeritas/mucf/interactor/DTMucfInteractor.hh index 39932258a8..ac6c5dec9f 100644 --- a/src/celeritas/mucf/interactor/DTMucfInteractor.hh +++ b/src/celeritas/mucf/interactor/DTMucfInteractor.hh @@ -59,6 +59,7 @@ class DTMucfInteractor Channel channel_{Channel::size_}; // Allocate space for secondary particles StackAllocator& allocate_; + // Number of secondaries per channel EnumArray num_secondaries_{ 3, // alpha_muon_neutron diff --git a/src/celeritas/mucf/model/detail/MucfMaterialInserter.cc b/src/celeritas/mucf/model/detail/MucfMaterialInserter.cc index f4d9930249..e6a7eaac1b 100644 --- a/src/celeritas/mucf/model/detail/MucfMaterialInserter.cc +++ b/src/celeritas/mucf/model/detail/MucfMaterialInserter.cc @@ -20,12 +20,12 @@ MucfMaterialInserter::MucfMaterialInserter(HostVal* host_data, inp::MucfPhysics const& data) : mucfmatid_to_matid_(&host_data->mucfmatid_to_matid) , isotopic_fractions_(&host_data->isotopic_fractions) - , cycle_times_(&host_data->cycle_times) + , cycle_rates_(&host_data->cycle_rates) , data_(data) { CELER_EXPECT(data_); - // Initialize interpolators for cycle time tables + // Initialize interpolators for cycle rate tables for (auto const& cycle_data : data_.cycle_rates) { // Use emplace to avoid copy/move of InterpolatorHelper objects @@ -43,7 +43,7 @@ bool MucfMaterialInserter::operator()(MaterialView const& material) using LhdArray = EquilibrateDensitiesSolver::LhdArray; MaterialFractionsArray isotopic_fractions; - CycleTimesArray cycle_times; + CycleRatesArray cycle_rates; LhdArray lhd_densities{}; auto from_mass_number = [&](AtomicMassNumber mass) -> MucfIsotope { @@ -94,19 +94,19 @@ bool MucfMaterialInserter::operator()(MaterialView const& material) if (lhd_densities[MucfIsotope::deuterium]) { - cycle_times[MucfMuonicMolecule::deuterium_deuterium] + cycle_rates[MucfMuonicMolecule::deuterium_deuterium] = this->calc_dd_cycle(equilibrium_densities, material.temperature()); } if (lhd_densities[MucfIsotope::tritium]) { - cycle_times[MucfMuonicMolecule::tritium_tritium] = this->calc_tt_cycle( + cycle_rates[MucfMuonicMolecule::tritium_tritium] = this->calc_tt_cycle( equilibrium_densities, material.temperature()); } if (lhd_densities[MucfIsotope::deuterium] && lhd_densities[MucfIsotope::tritium]) { - cycle_times[MucfMuonicMolecule::deuterium_tritium] + cycle_rates[MucfMuonicMolecule::deuterium_tritium] = this->calc_dt_cycle(equilibrium_densities, material.temperature()); } @@ -114,7 +114,7 @@ bool MucfMaterialInserter::operator()(MaterialView const& material) // Add muCF material to the model's host/device data mucfmatid_to_matid_.push_back(material.material_id()); isotopic_fractions_.push_back(std::move(isotopic_fractions)); - cycle_times_.push_back(std::move(cycle_times)); + cycle_rates_.push_back(std::move(cycle_rates)); //! \todo Store mean atom spin flip and transfer times @@ -123,7 +123,7 @@ bool MucfMaterialInserter::operator()(MaterialView const& material) //---------------------------------------------------------------------------// /*! - * Calculate dd muonic molecules cycle times. + * Calculate dd muonic molecules cycle rates. * * F = 1/2 and F = 3/2 are the reactive spin states for dd fusion. */ @@ -143,10 +143,8 @@ MucfMaterialInserter::calc_dd_cycle(EquilibriumArray const& eq_dens, = this->interpolator(CTT::deuterium_deuterium, HalfSpinInt{3}); MoleculeCycles result; - result[0] = real_type{1} - / (dd_dens * dd_1_over_2_interpolate(temperature)); // F = 1/2 - result[1] = real_type{1} - / (dd_dens * dd_3_over_2_interpolate(temperature)); // F = 3/2 + result[0] = (dd_dens * dd_1_over_2_interpolate(temperature)); // F = 1/2 + result[1] = (dd_dens * dd_3_over_2_interpolate(temperature)); // F = 3/2 CELER_ENSURE(result[0] >= 0 && result[1] >= 0); return result; @@ -154,7 +152,7 @@ MucfMaterialInserter::calc_dd_cycle(EquilibriumArray const& eq_dens, //---------------------------------------------------------------------------// /*! - * Calculate dt muonic molecules cycle times. + * Calculate dt muonic molecules cycle rates. * * F = 0 and F = 1 are the reactive spin states for dt fusion. */ @@ -189,14 +187,12 @@ MucfMaterialInserter::calc_dt_cycle(EquilibriumArray const& eq_dens, // Interpolate over rates, store final cycle time (1/rate) MoleculeCycles result; - result[0] = real_type{1} - / (hd_dens * hd0_interpolate(temperature) - + dd_dens * dd0_interpolate(temperature) - + dt_dens * dt0_interpolate(temperature)); // F = 0 - result[1] = real_type{1} - / (hd_dens * hd1_interpolate(temperature) - + dd_dens * dd1_interpolate(temperature) - + dt_dens * dt1_interpolate(temperature)); // F = 1 + result[0] = hd_dens * hd0_interpolate(temperature) + + dd_dens * dd0_interpolate(temperature) + + dt_dens * dt0_interpolate(temperature); // F = 0 + result[1] = hd_dens * hd1_interpolate(temperature) + + dd_dens * dd1_interpolate(temperature) + + dt_dens * dt1_interpolate(temperature); // F = 1 CELER_ENSURE(result[0] >= 0 && result[1] >= 0); return result; @@ -204,7 +200,7 @@ MucfMaterialInserter::calc_dt_cycle(EquilibriumArray const& eq_dens, //---------------------------------------------------------------------------// /*! - * Calculate tt muonic molecules cycle times. + * Calculate tt muonic molecules cycle rates. * * F = 1/2 is the only reactive spin state for tt fusion. */ @@ -221,7 +217,7 @@ MucfMaterialInserter::calc_tt_cycle(EquilibriumArray const& eq_dens, = this->interpolator(CTT::tritium_tritium, HalfSpinInt{1}); MoleculeCycles result; - result[0] = real_type{1} / (tt_dens * tt_interpolate(temperature)); + result[0] = tt_dens * tt_interpolate(temperature); // F = 1/2 CELER_ENSURE(result[0] >= 0 && result[1] == 0); return result; diff --git a/src/celeritas/mucf/model/detail/MucfMaterialInserter.hh b/src/celeritas/mucf/model/detail/MucfMaterialInserter.hh index 113605bb57..f5ed5cb1ec 100644 --- a/src/celeritas/mucf/model/detail/MucfMaterialInserter.hh +++ b/src/celeritas/mucf/model/detail/MucfMaterialInserter.hh @@ -47,7 +47,7 @@ class MucfMaterialInserter private: using MoleculeCycles = Array; - using CycleTimesArray = EnumArray; + using CycleRatesArray = EnumArray; using EquilibriumArray = EquilibrateDensitiesSolver::EquilibriumArray; using MaterialFractionsArray = EnumArray; using AtomicMassNumber = AtomicNumber; @@ -58,10 +58,10 @@ class MucfMaterialInserter //// DATA //// // DTMixMucfModel host data references populated by operator() - CollectionBuilder mucfmatid_to_matid_; - CollectionBuilder + CollectionBuilder mucfmatid_to_matid_; + CollectionBuilder isotopic_fractions_; - CollectionBuilder cycle_times_; + CollectionBuilder cycle_rates_; // Const data std::map const mass_isotope_map_{ {AtomicMassNumber{1}, MucfIsotope::protium}, @@ -73,15 +73,15 @@ class MucfMaterialInserter //// HELPER FUNCTIONS //// - // Calculate mean fusion cycle times for dd muonic molecules + // Calculate mean fusion cycle rates for dd muonic molecules Array calc_dd_cycle(EquilibriumArray const& eq_dens, real_type const temperature); - // Calculate mean fusion cycle times for dt muonic molecules + // Calculate mean fusion cycle rates for dt muonic molecules Array calc_dt_cycle(EquilibriumArray const& eq_dens, real_type const temperature); - // Calculate mean fusion cycle times for tt muonic molecules + // Calculate mean fusion cycle rates for tt muonic molecules Array calc_tt_cycle(EquilibriumArray const& eq_dens, real_type const temperature); diff --git a/test/celeritas/CMakeLists.txt b/test/celeritas/CMakeLists.txt index 1fe4465fb0..dbb13fe945 100644 --- a/test/celeritas/CMakeLists.txt +++ b/test/celeritas/CMakeLists.txt @@ -151,10 +151,12 @@ celeritas_add_test(decay/MuDecay.test.cc ${_needs_double}) #-----------------------------------------------------------------------------# # MuCF +celeritas_add_test(mucf/ChannelSelector.test.cc) celeritas_add_test(mucf/DDMucfInteractor.test.cc) celeritas_add_test(mucf/DTMucfInteractor.test.cc) celeritas_add_test(mucf/DTMixMucfModel.test.cc) celeritas_add_test(mucf/MuonicAtomSelector.test.cc) +celeritas_add_test(mucf/MuonicMoleculeSelector.test.cc) #-----------------------------------------------------------------------------# # EM diff --git a/test/celeritas/mucf/ChannelSelector.test.cc b/test/celeritas/mucf/ChannelSelector.test.cc new file mode 100644 index 0000000000..442c79aadf --- /dev/null +++ b/test/celeritas/mucf/ChannelSelector.test.cc @@ -0,0 +1,356 @@ +//------------------------------- -*- C++ -*- -------------------------------// +// Copyright Celeritas contributors: see top-level COPYRIGHT file for details +// SPDX-License-Identifier: (Apache-2.0 OR MIT) +//---------------------------------------------------------------------------// +//! \file celeritas/mucf/ChannelSelector.test.cc +//---------------------------------------------------------------------------// +#include "corecel/cont/Range.hh" +#include "corecel/random/DiagnosticRngEngine.hh" +#include "celeritas/mucf/executor/detail/DDChannelSelector.hh" +#include "celeritas/mucf/executor/detail/DTChannelSelector.hh" +#include "celeritas/mucf/executor/detail/TTChannelSelector.hh" + +#include "celeritas_test.hh" + +namespace celeritas +{ +namespace detail +{ +namespace test +{ +using namespace ::celeritas::test; +//---------------------------------------------------------------------------// +// TEST HARNESS +//---------------------------------------------------------------------------// + +class ChannelSelectorTest : public Test +{ + protected: + using Engine = DiagnosticRngEngine; + + void SetUp() override { rng_.reset_count(); } + + Engine& rng() { return rng_; } + + // Sticking fraction between the two dd --> 3He channels + double dd_sticking_fraction() { return 0.122; } + // Sticking fraction for dt + double dt_sticking_fraction() { return 0.00857; } + // Sticking fraction for tt + double tt_sticking_fraction() { return 0.14; } + + // Calculate dd --> 3He channel probability from the branching ratio + double he3_probability(double branching_ratio) + { + return branching_ratio / (branching_ratio + 1); + } + + // Calculate sigma for the statistical tests + double calc_sigma(double num_samples, double success_prob) + { + return std::sqrt(num_samples * success_prob * (1 - success_prob)); + } + + private: + Engine rng_; +}; + +//---------------------------------------------------------------------------// +// TESTS +//---------------------------------------------------------------------------// + +TEST_F(ChannelSelectorTest, dd_channel_low_temperature) +{ + // dd fusion at T < 50 K: branching_ratio = 1 + double const temperature = 30.0; + double const branching_ratio = 1.0; + double const he3_probability = this->he3_probability(branching_ratio); + double const sticking_fraction = this->dd_sticking_fraction(); + + DDChannelSelector select_channel(temperature); + + int num_samples = 100000; + int helium3_count = 0; + int muonichelium3_count = 0; + int tritium_count = 0; + + for ([[maybe_unused]] auto i : range(num_samples)) + { + auto channel = select_channel(this->rng()); + if (channel == DDMucfInteractor::Channel::helium3_muon_neutron) + { + helium3_count++; + } + else if (channel == DDMucfInteractor::Channel::muonichelium3_neutron) + { + muonichelium3_count++; + } + else if (channel == DDMucfInteractor::Channel::tritium_muon_proton) + { + tritium_count++; + } + else + { + FAIL() << "Unexpected channel selected"; + } + } + + EXPECT_EQ(num_samples, helium3_count + muonichelium3_count + tritium_count); + + if constexpr (CELERITAS_REAL_TYPE == CELERITAS_REAL_TYPE_DOUBLE) + { + double const num_samples_d = static_cast(num_samples); + + double const expected_tritium_count = num_samples_d + * (1 - he3_probability); + double const expected_helium3_count = num_samples_d * he3_probability + * (1 - sticking_fraction); + double const expected_muonichelium3_count + = num_samples_d * he3_probability * sticking_fraction; + double const tolerance + = 3 * this->calc_sigma(num_samples_d, he3_probability); + + EXPECT_NEAR(expected_tritium_count, tritium_count, tolerance); + EXPECT_NEAR(expected_helium3_count, helium3_count, tolerance); + EXPECT_NEAR( + expected_muonichelium3_count, muonichelium3_count, tolerance); + } +} + +//---------------------------------------------------------------------------// +TEST_F(ChannelSelectorTest, dd_channel_mid_temperature) +{ + // DD fusion at 50 < T < 100 K: branching_ratio = 1.0088 * (T - 50) = 25.22 + double const temperature = 75.0; + double const branching_ratio = 1.0088 * (temperature - 50); + double const he3_probability = this->he3_probability(branching_ratio); + + DDChannelSelector select_channel(temperature); + + size_type const num_samples = 10000; + size_type he3_total_count = 0; + size_type tritium_count = 0; + + for ([[maybe_unused]] auto i : range(num_samples)) + { + auto channel = select_channel(this->rng()); + if (channel == DDMucfInteractor::Channel::helium3_muon_neutron + || channel == DDMucfInteractor::Channel::muonichelium3_neutron) + { + he3_total_count++; + } + else if (channel == DDMucfInteractor::Channel::tritium_muon_proton) + { + tritium_count++; + } + } + + EXPECT_EQ(num_samples, he3_total_count + tritium_count); + + if constexpr (CELERITAS_REAL_TYPE == CELERITAS_REAL_TYPE_DOUBLE) + { + double const num_samples_d = static_cast(num_samples); + double const he3_total_count_d = static_cast(he3_total_count); + + double const expected_he3_count = num_samples_d * he3_probability; + double const tolerance + = 3 * this->calc_sigma(num_samples_d, he3_probability); + + EXPECT_NEAR(expected_he3_count, he3_total_count_d, tolerance); + } +} + +//---------------------------------------------------------------------------// +TEST_F(ChannelSelectorTest, dd_channel_high_temperature) +{ + // DD fusion at T < 300 K: branching_ratio = 1.44 + double const temperature = 300; + double const branching_ratio = 1.44; + double const he3_probability = this->he3_probability(branching_ratio); + + DDChannelSelector select_channel(temperature); + + size_type const num_samples = 10000; + size_type he3_total_count = 0; + size_type tritium_count = 0; + + for ([[maybe_unused]] auto i : range(num_samples)) + { + auto channel = select_channel(this->rng()); + if (channel == DDMucfInteractor::Channel::helium3_muon_neutron + || channel == DDMucfInteractor::Channel::muonichelium3_neutron) + { + he3_total_count++; + } + else if (channel == DDMucfInteractor::Channel::tritium_muon_proton) + { + tritium_count++; + } + } + + EXPECT_EQ(num_samples, he3_total_count + tritium_count); + + if constexpr (CELERITAS_REAL_TYPE == CELERITAS_REAL_TYPE_DOUBLE) + { + double const num_samples_d = static_cast(num_samples); + double const he3_total_count_d = static_cast(he3_total_count); + + double const expected_he3_count = num_samples_d * he3_probability; + // 3 sigma tolerance + double const tolerance + = 3 * this->calc_sigma(num_samples_d, he3_probability); + + EXPECT_NEAR(expected_he3_count, he3_total_count_d, tolerance); + } +} + +//---------------------------------------------------------------------------// +TEST_F(ChannelSelectorTest, dd_sticking_fraction_within_he3) +{ + // Test that when He3 channel is selected, sticking fraction is 12.2% + double const temperature = 300; + double const sticking_fraction = this->dd_sticking_fraction(); + + DDChannelSelector select_channel(temperature); + + size_type const num_samples = 10000; + size_type helium3_count = 0; + size_type muonichelium3_count = 0; + + // Only count He3 channels + for ([[maybe_unused]] auto i : range(num_samples)) + { + auto channel = select_channel(this->rng()); + if (channel == DDMucfInteractor::Channel::helium3_muon_neutron) + { + helium3_count++; + } + else if (channel == DDMucfInteractor::Channel::muonichelium3_neutron) + { + muonichelium3_count++; + } + } + + size_type total_he3 = helium3_count + muonichelium3_count; + EXPECT_GT(total_he3, 0); + + if constexpr (CELERITAS_REAL_TYPE == CELERITAS_REAL_TYPE_DOUBLE) + { + double const total_he3_d = static_cast(total_he3); + double const muonichelium3_count_d + = static_cast(muonichelium3_count); + + double const expected_muonichelium3 = total_he3_d * sticking_fraction; + // 3 sigma tolerance + double const tolerance + = 3 * this->calc_sigma(total_he3_d, sticking_fraction); + + EXPECT_NEAR(expected_muonichelium3, muonichelium3_count_d, tolerance); + } +} + +//---------------------------------------------------------------------------// +TEST_F(ChannelSelectorTest, dt_channel) +{ + // DT fusion: ~0.8% sticking + DTChannelSelector select_channel; + + size_type const num_samples = 100000; + size_type alpha_count = 0; + size_type muonicalpha_count = 0; + + for ([[maybe_unused]] auto i : range(num_samples)) + { + auto channel = select_channel(this->rng()); + if (channel == DTMucfInteractor::Channel::alpha_muon_neutron) + { + alpha_count++; + } + else if (channel == DTMucfInteractor::Channel::muonicalpha_neutron) + { + muonicalpha_count++; + } + else + { + FAIL() << "Unexpected channel selected"; + } + } + + EXPECT_EQ(num_samples, alpha_count + muonicalpha_count); + + if constexpr (CELERITAS_REAL_TYPE == CELERITAS_REAL_TYPE_DOUBLE) + { + double const num_samples_d = static_cast(num_samples); + double const alpha_count_d = static_cast(alpha_count); + double const muonicalpha_count_d + = static_cast(muonicalpha_count); + + double const sticking_fraction = this->dt_sticking_fraction(); + double const expected_muonicalpha_count = num_samples_d + * sticking_fraction; + double const expected_alpha_count = num_samples_d + * (1 - sticking_fraction); + // 3 sigma tolerance + double const tolerance + = 3 * this->calc_sigma(num_samples_d, sticking_fraction); + + EXPECT_NEAR(expected_muonicalpha_count, muonicalpha_count_d, tolerance); + EXPECT_NEAR(expected_alpha_count, alpha_count_d, tolerance); + } +} + +//---------------------------------------------------------------------------// +TEST_F(ChannelSelectorTest, tt_channel) +{ + // TT fusion: ~14% sticking + TTChannelSelector select_channel; + + size_type const num_samples = 10000; + size_type alpha_count = 0; + size_type muonicalpha_count = 0; + + for ([[maybe_unused]] auto i : range(num_samples)) + { + auto channel = select_channel(this->rng()); + if (channel == TTMucfInteractor::Channel::alpha_muon_neutron_neutron) + { + alpha_count++; + } + else if (channel + == TTMucfInteractor::Channel::muonicalpha_neutron_neutron) + { + muonicalpha_count++; + } + else + { + FAIL() << "Unexpected channel selected"; + } + } + + EXPECT_EQ(num_samples, alpha_count + muonicalpha_count); + + if constexpr (CELERITAS_REAL_TYPE == CELERITAS_REAL_TYPE_DOUBLE) + { + double const num_samples_d = static_cast(num_samples); + double const alpha_count_d = static_cast(alpha_count); + double const muonicalpha_count_d + = static_cast(muonicalpha_count); + + double const sticking_fraction = this->tt_sticking_fraction(); + double const expected_muonicalpha_count = num_samples_d + * sticking_fraction; + double const expected_alpha_count = num_samples_d + * (1 - sticking_fraction); + + double const tolerance + = 3 * this->calc_sigma(num_samples_d, sticking_fraction); + + EXPECT_NEAR(expected_muonicalpha_count, muonicalpha_count_d, tolerance); + EXPECT_NEAR(expected_alpha_count, alpha_count_d, tolerance); + } +} + +//---------------------------------------------------------------------------// +} // namespace test +} // namespace detail +} // namespace celeritas diff --git a/test/celeritas/mucf/DTMixMucfModel.test.cc b/test/celeritas/mucf/DTMixMucfModel.test.cc index f51f74f28c..8701727bef 100644 --- a/test/celeritas/mucf/DTMixMucfModel.test.cc +++ b/test/celeritas/mucf/DTMixMucfModel.test.cc @@ -69,9 +69,9 @@ TEST_F(DTMixMucfModelTest, data) // Check sotopic fractions // Single material with 50/50 d and t fractions EXPECT_SOFT_EQ( - 0.5, data.isotopic_fractions[MuCfMatId{0}][MucfIsotope::deuterium]); + 0.5, data.isotopic_fractions[MucfMatId{0}][MucfIsotope::deuterium]); EXPECT_SOFT_EQ( - 0.5, data.isotopic_fractions[MuCfMatId{0}][MucfIsotope::tritium]); + 0.5, data.isotopic_fractions[MucfMatId{0}][MucfIsotope::tritium]); // Cycle times are in seconds // DD (reactivity of F = 3/2 is almost negligible, with huge cycle times) @@ -83,22 +83,23 @@ TEST_F(DTMixMucfModelTest, data) // TT static double const expected_tt_1_over_2_cycle_time{1.4056833511329384e-06}; - auto const& cycles = data.cycle_times; + auto const& rates = data.cycle_rates; // DD cycle times EXPECT_SOFT_EQ(expected_dd_1_over_2_cycle_time, - cycles[MuCfMatId{0}][Molecule::deuterium_deuterium][0]); + 1 / rates[MucfMatId{0}][Molecule::deuterium_deuterium][0]); EXPECT_SOFT_EQ(expected_dd_3_over_2_cycle_time, - cycles[MuCfMatId{0}][Molecule::deuterium_deuterium][1]); + 1 / rates[MucfMatId{0}][Molecule::deuterium_deuterium][1]); // DT cycle times EXPECT_SOFT_EQ(expected_dt_0_cycle_time, - cycles[MuCfMatId{0}][Molecule::deuterium_tritium][0]); + 1 / rates[MucfMatId{0}][Molecule::deuterium_tritium][0]); EXPECT_SOFT_EQ(expected_dt_1_cycle_time, - cycles[MuCfMatId{0}][Molecule::deuterium_tritium][1]); + 1 / rates[MucfMatId{0}][Molecule::deuterium_tritium][1]); // TT cycle times EXPECT_SOFT_EQ(expected_tt_1_over_2_cycle_time, - cycles[MuCfMatId{0}][Molecule::tritium_tritium][0]); - EXPECT_SOFT_EQ(0, cycles[MuCfMatId{0}][Molecule::tritium_tritium][1]); + 1 / rates[MucfMatId{0}][Molecule::tritium_tritium][0]); + EXPECT_SOFT_EQ(numeric_limits::infinity(), + 1 / rates[MucfMatId{0}][Molecule::tritium_tritium][1]); } //---------------------------------------------------------------------------// diff --git a/test/celeritas/mucf/MuonicMoleculeSelector.test.cc b/test/celeritas/mucf/MuonicMoleculeSelector.test.cc new file mode 100644 index 0000000000..256b872ffa --- /dev/null +++ b/test/celeritas/mucf/MuonicMoleculeSelector.test.cc @@ -0,0 +1,161 @@ +//------------------------------- -*- C++ -*- -------------------------------// +// Copyright Celeritas contributors: see top-level COPYRIGHT file for details +// SPDX-License-Identifier: (Apache-2.0 OR MIT) +//---------------------------------------------------------------------------// +//! \file celeritas/mucf/MuonicMoleculeSelector.test.cc +//---------------------------------------------------------------------------// +#include "celeritas/mucf/executor/detail/MuonicMoleculeSelector.hh" + +#include "corecel/cont/Range.hh" +#include "corecel/io/Logger.hh" +#include "corecel/random/DiagnosticRngEngine.hh" +#include "celeritas/mucf/Types.hh" +#include "celeritas/mucf/executor/detail/MuonicAtomSelector.hh" +#include "celeritas/mucf/executor/detail/MuonicAtomSpinSelector.hh" + +#include "MucfInteractorHostTestBase.hh" +#include "celeritas_test.hh" + +namespace celeritas +{ +namespace detail +{ +namespace test +{ +using namespace ::celeritas::test; +//---------------------------------------------------------------------------// +// TEST HARNESS +//---------------------------------------------------------------------------// + +class MuonicMoleculeSelectorTest + : public ::celeritas::test::MucfInteractorHostTestBase +{ + protected: + using Base = ::celeritas::test::MucfInteractorHostTestBase; + using Engine = DiagnosticRngEngine; + using MMM = MucfMuonicMolecule; + using MMA = MucfMuonicAtom; + using HalfSpinInt = MuonicMoleculeSelector::HalfSpinInt; + using CycleRatesArray = MuonicMoleculeSelector::CycleRatesArray; + + void SetUp() override + { + data_ = this->host_data(); + rng_.reset_count(); + + mucfmatid_ = MucfMatId{0}; + cycle_rates_ = data_.cycle_rates[mucfmatid_]; + } + + protected: + HostCRef data_; + Engine rng_; + MucfMatId mucfmatid_; + CycleRatesArray cycle_rates_; +}; + +//---------------------------------------------------------------------------// +// TESTS +//---------------------------------------------------------------------------// +TEST_F(MuonicMoleculeSelectorTest, TEST_IF_CELERITAS_DOUBLE(dd_spin_selection)) +{ + CycleRatesArray rates; + + // Make DT and TT rates zero + rates[MMM::deuterium_tritium] = {0, 0}; + rates[MMM::tritium_tritium] = {0, 0}; + + // DD, F = 1/2 + rates[MMM::deuterium_deuterium] = {1e6, 0}; + auto const dd_1_result + = MuonicMoleculeSelector(MMA::deuterium, HalfSpinInt{1}, rates)(rng_); + + EXPECT_EQ(dd_1_result.molecule, MMM::deuterium_deuterium); + EXPECT_SOFT_EQ(dd_1_result.cycle_time, 1.9989533630516e-06); + + // DD, F = 3/2 + rates[MMM::deuterium_deuterium] = {0, 1e6}; + auto const dd_3_result + = MuonicMoleculeSelector(MMA::deuterium, HalfSpinInt{3}, rates)(rng_); + + EXPECT_EQ(dd_3_result.molecule, MMM::deuterium_deuterium); + EXPECT_SOFT_EQ(dd_3_result.cycle_time, 1.8031326676554e-07); +} + +//---------------------------------------------------------------------------// +TEST_F(MuonicMoleculeSelectorTest, TEST_IF_CELERITAS_DOUBLE(dt_spin_selection)) +{ + CycleRatesArray rates; + + // Make DD and TT rates zero + rates[MMM::deuterium_deuterium] = {0, 0}; + rates[MMM::tritium_tritium] = {0, 0}; + + // DT, F = 0 + rates[MMM::deuterium_tritium] = {1e6, 0}; + auto const dt_0_result + = MuonicMoleculeSelector(MMA::tritium, HalfSpinInt{0}, rates)(rng_); + + EXPECT_EQ(dt_0_result.molecule, MMM::deuterium_tritium); + EXPECT_SOFT_EQ(dt_0_result.cycle_time, 1.9989533630516e-06); + + // DT, F = 1 + rates[MMM::deuterium_tritium] = {0, 1e6}; + auto const dt_2_result + = MuonicMoleculeSelector(MMA::tritium, HalfSpinInt{2}, rates)(rng_); + + EXPECT_EQ(dt_2_result.molecule, MMM::deuterium_tritium); + EXPECT_SOFT_EQ(dt_2_result.cycle_time, 1.8031326676554e-07); +} + +//---------------------------------------------------------------------------// +TEST_F(MuonicMoleculeSelectorTest, TEST_IF_CELERITAS_DOUBLE(tt_spin_selection)) +{ + CycleRatesArray rates; + + // TT, F = 1/2 + rates[MMM::deuterium_deuterium] = {0, 0}; + rates[MMM::deuterium_tritium] = {0, 0}; + rates[MMM::tritium_tritium] = {1e6, 0}; + + auto const tt_1_result + = MuonicMoleculeSelector(MMA::tritium, HalfSpinInt{1}, rates)(rng_); + EXPECT_EQ(tt_1_result.molecule, MMM::tritium_tritium); + EXPECT_SOFT_EQ(tt_1_result.cycle_time, 1.9989533630516e-06); +} + +//---------------------------------------------------------------------------// +TEST_F(MuonicMoleculeSelectorTest, model_data) +{ + using MoleculeCountArray = EnumArray; + + // Mimic executor behavior: form atom, then form molecule + size_type const num_samples = 10000; + auto const& rates = data_.cycle_rates[mucfmatid_]; + + MoleculeCountArray molecule_counts; + for ([[maybe_unused]] auto i : range(num_samples)) + { + // Form muonic atom + auto muonic_atom = MuonicAtomSelector( + data_.isotopic_fractions[mucfmatid_][MucfIsotope::deuterium])(rng_); + auto atom_spin = MuonicAtomSpinSelector(muonic_atom)(rng_); + + // Form molecule that will call determine the interactor + auto result + = MuonicMoleculeSelector(muonic_atom, atom_spin, rates)(rng_); + molecule_counts[result.molecule]++; + } + + // From Acceleron's simulation with equivalent material: 50/50 d/t at 300K, + // the fraction of each fusion call is: DD: 20.1%, DT: 78.9%, TT: 0.9%. + // The counts on this test (20%, 79%, 0.2%) are close to Acceleron's. + EXPECT_EQ(2044, molecule_counts[MMM::deuterium_deuterium]); + EXPECT_EQ(7940, molecule_counts[MMM::deuterium_tritium]); + EXPECT_EQ(16, molecule_counts[MMM::tritium_tritium]); +} + +//---------------------------------------------------------------------------// +} // namespace test +} // namespace detail +} // namespace celeritas