diff --git a/doc/CMakeLists.txt b/doc/CMakeLists.txt index eb2d012e88..ad297309b8 100644 --- a/doc/CMakeLists.txt +++ b/doc/CMakeLists.txt @@ -263,6 +263,12 @@ file(GLOB _DOXYGEN_SOURCE "${PROJECT_SOURCE_DIR}/src/celeritas/em/msc/detail/*.hh" "${PROJECT_SOURCE_DIR}/src/celeritas/em/process/*.hh" "${PROJECT_SOURCE_DIR}/src/celeritas/em/xs/*.hh" + "${PROJECT_SOURCE_DIR}/src/celeritas/mucf/data/*.hh" + "${PROJECT_SOURCE_DIR}/src/celeritas/mucf/executor/*.hh" + "${PROJECT_SOURCE_DIR}/src/celeritas/mucf/interactor/*.hh" + "${PROJECT_SOURCE_DIR}/src/celeritas/mucf/model/*.hh" + "${PROJECT_SOURCE_DIR}/src/celeritas/mucf/model/detail/*.hh" + "${PROJECT_SOURCE_DIR}/src/celeritas/mucf/process/*.hh" "${PROJECT_SOURCE_DIR}/src/celeritas/neutron/interactor/*.hh" "${PROJECT_SOURCE_DIR}/src/celeritas/neutron/interactor/detail/*.hh" "${PROJECT_SOURCE_DIR}/src/celeritas/neutron/model/*.hh" diff --git a/doc/implementation/mucf-physics.rst b/doc/implementation/mucf-physics.rst index ec6d81f84e..75e13de3bf 100644 --- a/doc/implementation/mucf-physics.rst +++ b/doc/implementation/mucf-physics.rst @@ -72,8 +72,8 @@ sticking factor and the fusion cycle time are the main conditions that define how many fusion cycles a muon can undergo. The fusion cycle time depends on the d-t mixture, its temperature, and on the final spin of the molecule. Only muonic molecules where the total spin :math:`F = I_N \pm 1/2` is on, or has a -projection onto the total angular momentum J = 1 are reactive. The spin states -of the three possible muonic molecules are summarized in table +projection onto the total angular momentum :math:`J = 1` are reactive. The spin +states of the three possible muonic molecules are summarized in table :numref:`muon_spin_states`. @@ -108,9 +108,36 @@ enabling the ``mucf_physics`` option in Geant4 integration ------------------ -For integration interfaces, if ``mucf_physics`` option in -:cpp:class:`celeritas::ext::GeantPhysicsOptions` is enabled, the muon-catalyzed -fusion data is constructed when the ``G4MuonMinusAtomicCapture`` process is -registered. +For integration interfaces, enabling the ``mucf_physics`` option in +:cpp:class:`celeritas::ext::GeantPhysicsOptions` will check if the +``G4MuonMinusAtomicCapture`` process is registered in the Geant4's Physics List. +If the process is present, the :cpp:class:`celeritas::inp::MucfPhysics` will be +populated, and the :cpp:class:`celeritas::MucfProcess` will be initialized. -.. todo:: Add process/model/executor details +Code implementation +=================== + +The :cpp:class:`celeritas::MucfProcess` process has only the +:cpp:class:`celeritas::DTMixMucfModel` attached to it, responsible for +deuterium-tritium mixtures. It can simulate materials from near absolute zero to +1500 kelvin. It is an *at rest* model that encompasses the full cycle---atom +formation, molecule formation, and fusion. + +.. doxygenclass:: celeritas::MucfProcess +.. doxygenclass:: celeritas::DTMixMucfModel + +Most of the data is material-dependent, being calculated and cached during model +construction. All of the cached quantities are calculated and added to +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`, with the Interactors reserved for +sampling final states of the outgoing secondaries. + +.. note:: Only reactive channels are implemented. + +.. doxygenclass:: celeritas::DDMucfInteractor +.. doxygenclass:: celeritas::DTMucfInteractor +.. doxygenclass:: celeritas::TTMucfInteractor diff --git a/src/celeritas/CMakeLists.txt b/src/celeritas/CMakeLists.txt index 4a108c8774..b0c122e52d 100644 --- a/src/celeritas/CMakeLists.txt +++ b/src/celeritas/CMakeLists.txt @@ -100,6 +100,8 @@ list(APPEND SOURCES mat/MaterialParams.cc mat/MaterialParamsOutput.cc mat/detail/Utils.cc + mucf/model/detail/MucfMaterialInserter.cc + mucf/process/MucfProcess.cc neutron/model/CascadeOptions.cc neutron/model/CascadeOptionsIO.json.cc neutron/process/NeutronElasticProcess.cc @@ -375,6 +377,7 @@ celeritas_polysource(em/model/CoulombScatteringModel) celeritas_polysource(geo/detail/BoundaryAction) celeritas_polysource(global/detail/KillActive) celeritas_polysource(global/detail/TrackSlotUtils) +celeritas_polysource(mucf/model/DTMixMucfModel) celeritas_polysource(neutron/model/ChipsNeutronElasticModel) celeritas_polysource(neutron/model/NeutronInelasticModel) celeritas_polysource(optical/model/AbsorptionModel) diff --git a/src/celeritas/Types.hh b/src/celeritas/Types.hh index c2bc81764a..c2cd0e9e63 100644 --- a/src/celeritas/Types.hh +++ b/src/celeritas/Types.hh @@ -227,25 +227,6 @@ enum class CylAxis size_ }; -//---------------------------------------------------------------------------// -//! Muon-catalyzed fusion atoms -enum class MucfMuonicAtom -{ - deuterium, - tritium, - size_ -}; - -//---------------------------------------------------------------------------// -//! Muon-catalyzed fusion molecules -enum class MucfMuonicMolecule -{ - deuterium_deuterium, - deuterium_tritium, - tritium_tritium, - size_ -}; - //---------------------------------------------------------------------------// // HELPER STRUCTS //---------------------------------------------------------------------------// diff --git a/src/celeritas/ext/RootInterfaceLinkDef.h b/src/celeritas/ext/RootInterfaceLinkDef.h index 3e1ac8f498..0425bc3110 100644 --- a/src/celeritas/ext/RootInterfaceLinkDef.h +++ b/src/celeritas/ext/RootInterfaceLinkDef.h @@ -51,10 +51,11 @@ #pragma link C++ class celeritas::inp::Grid+; #pragma link C++ class celeritas::inp::GridReflection+; #pragma link C++ class celeritas::inp::Interpolation+; -#pragma link C++ class celeritas::inp::MucfPhysics+; -#pragma link C++ class celeritas::inp::MucfCycleRate+; #pragma link C++ class celeritas::inp::MucfAtomTransferRate+; #pragma link C++ class celeritas::inp::MucfAtomSpinFlipRate+; +#pragma link C++ class celeritas::inp::MucfCycleRate+; +#pragma link C++ class celeritas::inp::MucfPhysics+; +#pragma link C++ class celeritas::inp::MucfScalars+; #pragma link C++ class celeritas::inp::MuPairProductionEnergyTransferTable+; #pragma link C++ class celeritas::inp::NoRoughness+; #pragma link C++ class celeritas::inp::OpticalPhysics+; @@ -82,6 +83,9 @@ // Quantities #pragma link C++ class celeritas::Quantity+; #pragma link C++ class celeritas::Quantity+; +#pragma link C++ class celeritas::Quantity+; +#pragma link C++ class celeritas::Quantity+; + // Event data used by Geant4/Celeritas offloading applications #pragma link C++ class celeritas::EventData+; diff --git a/src/celeritas/inp/MucfPhysics.cc b/src/celeritas/inp/MucfPhysics.cc index 9724a6e0b0..9a298eaa7b 100644 --- a/src/celeritas/inp/MucfPhysics.cc +++ b/src/celeritas/inp/MucfPhysics.cc @@ -10,6 +10,117 @@ namespace celeritas { namespace inp { +namespace +{ +//---------------------------------------------------------------------------// +/*! + * Muon energy CDF data for muon-catalyzed fusion. + */ +Grid mucf_muon_energy_cdf() +{ + Grid cdf; + cdf.interpolation.type = InterpolationType::cubic_spline; + + // Cumulative distribution data [unitless] + cdf.x = { + 0, + 0.04169381107491854, + 0.08664495114006499, + 0.14332247557003264, + 0.20456026058631915, + 0.2723127035830618, + 0.34136807817589576, + 0.41563517915309456, + 0.48990228013029324, + 0.5667752442996744, + 0.6306188925081434, + 0.6866449511400652, + 0.7309446254071662, + 0.7778501628664496, + 0.8104234527687297, + 0.8403908794788275, + 0.8618892508143323, + 0.8814332247557004, + 0.8970684039087949, + 0.903583061889251, + 1.0, + }; + + // Energy [keV] + cdf.y = { + 0, + 0.48850540675768084, + 0.8390389347819425, + 1.2521213482687141, + 1.7153033196164724, + 2.253638712180777, + 2.854653691809707, + 3.606073540073316, + 4.470346052913727, + 5.560291219507215, + 6.700556502915258, + 7.953772477101693, + 9.194596305637525, + 10.849180562221111, + 12.353474314071864, + 14.045888515617822, + 15.650634617544647, + 17.38079707555165, + 19.111008546659452, + 19.976130619913615, + 80.0, + }; + + CELER_ENSURE(cdf); + return cdf; +} + +//---------------------------------------------------------------------------// +/*! + * Cycle rate data for muon-catalyzed fusion. + */ +std::vector mucf_cycle_rates() +{ + //! \todo Implement + return {}; +} + +//---------------------------------------------------------------------------// +/*! + * Spin-flip data for muon-catalyzed fusion. + */ +std::vector mucf_atom_spin_flip_rates() +{ + //! \todo Implement + return {}; +} + +//---------------------------------------------------------------------------// +/*! + * Atom transfer data for muon-catalyzed fusion. + */ +std::vector mucf_atom_transfer_rates() +{ + //! \todo Implement + return {}; +} +//---------------------------------------------------------------------------// +} // namespace + +//---------------------------------------------------------------------------// +/*! + * Initialize with hardcoded values. + */ +MucfScalars MucfScalars::from_default() +{ + MucfScalars result; + result.protium = AmuMass{1.007825031898}; + result.deuterium = AmuMass{2.014101777844}; + result.tritium = AmuMass{3.016049281320}; + result.liquid_hydrogen_density = InvCcDensity{4.25e22}; + return result; +} + //---------------------------------------------------------------------------// /*! * Construct hardcoded muon-catalyzed fusion physics data. @@ -21,16 +132,14 @@ namespace inp MucfPhysics MucfPhysics::from_default() { MucfPhysics result; - - //! \todo Initialize hardcoded CDF data - //! \todo Initialize hardcoded cycle rate data - //! \todo Initialize hardcoded atom transfer data - //! \todo Initialize hardcoded spin flip data + result.scalars = MucfScalars::from_default(); + result.muon_energy_cdf = mucf_muon_energy_cdf(); + result.cycle_rates = mucf_cycle_rates(); + result.atom_transfer = mucf_atom_transfer_rates(); + result.atom_spin_flip = mucf_atom_spin_flip_rates(); // Temporary test dummy data to verify correct import { - result.muon_energy_cdf = Grid::from_constant(1.0); - MucfCycleRate dt_cycle; dt_cycle.molecule = MucfMuonicMolecule::deuterium_tritium; diff --git a/src/celeritas/inp/MucfPhysics.hh b/src/celeritas/inp/MucfPhysics.hh index 88c67a90e4..401e58b525 100644 --- a/src/celeritas/inp/MucfPhysics.hh +++ b/src/celeritas/inp/MucfPhysics.hh @@ -10,12 +10,43 @@ #include #include "corecel/inp/Grid.hh" -#include "celeritas/Types.hh" +#include "corecel/math/Quantity.hh" +#include "celeritas/UnitTypes.hh" +#include "celeritas/mucf/Types.hh" namespace celeritas { namespace inp { +//---------------------------------------------------------------------------// +/*! + * Muon-catalyzed fusion scalars. + * + * Default values are the same used by Acceleron. + */ +struct MucfScalars +{ + using AmuMass = Quantity; + using InvCcDensity = Quantity; + + // Atomic masses + AmuMass protium; //!< Protium atomic mass [AMU] + AmuMass deuterium; //!< Deuterium atomic mass [AMU] + AmuMass tritium; //!< Tritium atomic mass [AMU] + InvCcDensity liquid_hydrogen_density; //!< LHD unit [1/cm^3] + + //! Whether scalars have been defined + explicit operator bool() const + { + return protium > zero_quantity() && deuterium > zero_quantity() + && tritium > zero_quantity() + && liquid_hydrogen_density > zero_quantity(); + } + + // Initialize with hardcoded values + static MucfScalars from_default(); +}; + //---------------------------------------------------------------------------// /*! * Muon-catalyzed fusion mean cycle rate data. @@ -68,6 +99,9 @@ struct MucfCycleRate struct MucfAtomTransferRate { //! \todo Implement + + //! True if data is assigned \todo Implement + explicit operator bool() const { return true; } }; //---------------------------------------------------------------------------// @@ -91,6 +125,9 @@ struct MucfAtomTransferRate struct MucfAtomSpinFlipRate { //! \todo Implement + + //! True if data is assigned \todo Implement + explicit operator bool() const { return true; } }; //---------------------------------------------------------------------------// @@ -98,7 +135,8 @@ struct MucfAtomSpinFlipRate * Muon-catalyzed fusion physics options and data import. * * Minimum requirements for muon-catalyzed fusion: - * - Muon energy CDF data, required for sampling the outgoing muCF muon, and + * - Muon energy CDF data, required for sampling the outgoing muCF muon, + * and * - Mean cycle rate data for dd, dt, and tt muonic molecules. * * Muonic atom transfer and muonic atom spin flip are secondary effects and not @@ -109,15 +147,16 @@ struct MucfPhysics template using Vec = std::vector; - Grid muon_energy_cdf; //!< CDF for outgoing muCF muon + MucfScalars scalars; + Grid muon_energy_cdf; //!< CDF for sampling the outgoing muCF muon Vec cycle_rates; //!< Mean cycle rates for muonic molecules - Vec atom_transfer; //!< Muon atom transfer rates - Vec atom_spin_flip; //!< Muon atom spin flip rates + Vec atom_transfer; //!< Muonic atom transfer rates + Vec atom_spin_flip; //!< Muonic atom spin flip rates //! Whether muon-catalyzed fusion physics is enabled explicit operator bool() const { - return muon_energy_cdf && !cycle_rates.empty(); + return scalars && muon_energy_cdf && !cycle_rates.empty(); } //! Construct hardcoded muon-catalyzed fusion physics data diff --git a/src/celeritas/mucf/Types.hh b/src/celeritas/mucf/Types.hh new file mode 100644 index 0000000000..0a7c410dda --- /dev/null +++ b/src/celeritas/mucf/Types.hh @@ -0,0 +1,71 @@ +//------------------------------- -*- C++ -*- -------------------------------// +// Copyright Celeritas contributors: see top-level COPYRIGHT file for details +// SPDX-License-Identifier: (Apache-2.0 OR MIT) +//---------------------------------------------------------------------------// +//! \file celeritas/mucf/Types.hh +//---------------------------------------------------------------------------// +#pragma once + +#include "corecel/OpaqueId.hh" + +namespace celeritas +{ +//---------------------------------------------------------------------------// +// ENUMERATIONS +//---------------------------------------------------------------------------// +/*! + * Muonic atom selection from material data. This is *not* intended to be used + * by the transport loop. + */ +enum class MucfMuonicAtom +{ + deuterium, + tritium, + size_ +}; + +//---------------------------------------------------------------------------// +/*! + * Muonic molecule selection from material data. This is *not* intended to be + * used by the transport loop. + */ +enum class MucfMuonicMolecule +{ + deuterium_deuterium, + deuterium_tritium, + tritium_tritium, + size_ +}; + +//---------------------------------------------------------------------------// +/*! + * Enum for safely accessing hydrogen isoprotologues. + * + * Hydrogen isoprotologue molecules are: + * - Homonuclear: \f$ ^2H \f$, \f$ ^2d \f$, and \f$ ^2t \f$ + * - Heteronuclear: hd, ht, and dt. + * + * \note Muon-catalyzed fusion data is only applicable to a material with + * concentrations in thermodynamic equilibrium. This equilibrium is calculated + * at model construction from the material temperature and its h, d, and t + * fractions. + */ +enum class MucfIsoprotologueMolecule +{ + protium_protium, + protium_deuterium, + protium_tritium, + deuterium_tritium, + tritium_tritium, + size_ +}; + +//---------------------------------------------------------------------------// +// TYPE ALIASES +//---------------------------------------------------------------------------// + +//! Opaque index of a muCF material component +using MuCfMatId = OpaqueId; + +//---------------------------------------------------------------------------// +} // namespace celeritas diff --git a/src/celeritas/mucf/data/DTMixMucfData.hh b/src/celeritas/mucf/data/DTMixMucfData.hh new file mode 100644 index 0000000000..63affe3867 --- /dev/null +++ b/src/celeritas/mucf/data/DTMixMucfData.hh @@ -0,0 +1,113 @@ +//------------------------------- -*- C++ -*- -------------------------------// +// Copyright Celeritas contributors: see top-level COPYRIGHT file for details +// SPDX-License-Identifier: (Apache-2.0 OR MIT) +//---------------------------------------------------------------------------// +//! \file celeritas/mucf/data/DTMixMucfData.hh +//---------------------------------------------------------------------------// +#pragma once + +#include +#include + +#include "corecel/Macros.hh" +#include "corecel/cont/EnumArray.hh" +#include "corecel/grid/NonuniformGridData.hh" +#include "corecel/io/Join.hh" +#include "celeritas/Types.hh" +#include "celeritas/mucf/Types.hh" +#include "celeritas/phys/ParticleParams.hh" + +namespace celeritas +{ +//---------------------------------------------------------------------------// +/*! + * ParticleIds used by the \c DTMixMucfModel . + */ +struct MucfParticleIds +{ + //! Primary + ParticleId mu_minus; + + //!@{ + //! Elementary particles and nuclei + ParticleId neutron; + ParticleId proton; + ParticleId alpha; + ParticleId he3; + //!@} + + //!@{ + //! Muonic atoms + ParticleId muonic_hydrogen; + ParticleId muonic_deuteron; + ParticleId muonic_triton; + ParticleId muonic_alpha; + ParticleId muonic_he3; + //!@} + + //! Check whether all particles are assigned + CELER_FUNCTION explicit operator bool() const + { + return mu_minus && neutron && proton && alpha && he3 && muonic_hydrogen + && muonic_alpha && muonic_triton && muonic_he3; + } +}; + +//---------------------------------------------------------------------------// +/*! + * Data for for the \c DTMixMucfModel . + */ +template +struct DTMixMucfData +{ + template + using Items = Collection; + template + using MaterialItems = Collection; + using GridRecord = NonuniformGridRecord; + using CycleTimesArray = EnumArray>; + + //! Particle IDs + MucfParticleIds particles; + + //! Muon CDF energy grid for sampling outgoing muCF muons + GridRecord muon_energy_cdf; //! \todo Verify energy unit + Items reals; + + //!@{ + //! Material-dependent data calculated at model construction + //! \c PhysMatId indexed by \c MuCfMatId + MaterialItems mucfmatid_to_matid; + //! Cycle times per material: [mat_comp_id][muonic_molecule][spin_index] + MaterialItems cycle_times; //!< In [s] + //! \todo Add mean atom spin flip times + //! \todo Add mean atom transfer times + //!@} + + //! Check whether the data are assigned + explicit CELER_FUNCTION operator bool() const + { + return particles && muon_energy_cdf && !mucfmatid_to_matid.empty() + && !cycle_times.empty() + && (mucfmatid_to_matid.size() == cycle_times.size()); + } + + //! Assign from another set of data + template + DTMixMucfData& operator=(DTMixMucfData const& other) + { + CELER_EXPECT(other); + + //! \todo Finish implementation + this->particles = other.particles; + this->reals = other.reals; + this->muon_energy_cdf = other.muon_energy_cdf; + this->mucfmatid_to_matid = other.mucfmatid_to_matid; + this->cycle_times = other.cycle_times; + + return *this; + } +}; + +//---------------------------------------------------------------------------// +} // namespace celeritas diff --git a/src/celeritas/mucf/executor/DTMixMucfExecutor.hh b/src/celeritas/mucf/executor/DTMixMucfExecutor.hh new file mode 100644 index 0000000000..030d2c2bf9 --- /dev/null +++ b/src/celeritas/mucf/executor/DTMixMucfExecutor.hh @@ -0,0 +1,148 @@ +//------------------------------- -*- C++ -*- -------------------------------// +// Copyright Celeritas contributors: see top-level COPYRIGHT file for details +// SPDX-License-Identifier: (Apache-2.0 OR MIT) +//---------------------------------------------------------------------------// +//! \file celeritas/mucf/executor/DTMixMucfExecutor.hh +//---------------------------------------------------------------------------// +#pragma once + +#include "corecel/Assert.hh" +#include "celeritas/global/CoreTrackView.hh" +#include "celeritas/mat/ElementView.hh" +#include "celeritas/mucf/data/DTMixMucfData.hh" +#include "celeritas/mucf/interactor/DDMucfInteractor.hh" +#include "celeritas/mucf/interactor/DTMucfInteractor.hh" +#include "celeritas/mucf/interactor/TTMucfInteractor.hh" + +#include "detail/DDChannelSelector.hh" +#include "detail/DTChannelSelector.hh" +#include "detail/DTMixMuonicAtomSelector.hh" +#include "detail/DTMixMuonicMoleculeSelector.hh" +#include "detail/MuonicAtomSpinSelector.hh" +#include "detail/MuonicMoleculeSpinSelector.hh" +#include "detail/TTChannelSelector.hh" + +namespace celeritas +{ +//---------------------------------------------------------------------------// +struct DTMixMucfExecutor +{ + inline CELER_FUNCTION Interaction + operator()(celeritas::CoreTrackView const& track); + + NativeCRef data; +}; + +//---------------------------------------------------------------------------// +/*! + * Execute muon-catalyzed fusion for muonic dd, dt, or tt molecules. + */ +CELER_FUNCTION Interaction +DTMixMucfExecutor::operator()(celeritas::CoreTrackView const& track) +{ + auto phys_step_view = track.physics_step(); + auto elcomp_id = phys_step_view.element(); + CELER_ASSERT(elcomp_id); + + auto element = track.material().material_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 + + // Form d or t muonic atom + detail::DTMixMuonicAtomSelector form_atom; + auto muonic_atom = form_atom(rng); + + // Select atom spin via a helper class + detail::MuonicAtomSpinSelector select_atom_spin(muonic_atom); + + // { + // Competing at-rest processes which add to the total track time + //! \todo Muonic atom transfer + //! \todo Muonic atom spin flip + // } + + // 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); + + // Find muCF material ID from PhysMatId + // Make this a View if ever used beyond this executor + auto find = [&](PhysMatId matid) -> MuCfMatId { + CELER_EXPECT(matid); + for (auto i : range(data.mucfmatid_to_matid.size())) + { + if (auto const comp_id = MuCfMatId{i}; + data.mucfmatid_to_matid[comp_id] == matid) + { + return comp_id; + } + } + // MuCF material ID not found + return MuCfMatId{}; + }; + + // Load cycle time for the selected molecule + auto const mucf_matid = find(track.material().material_id()); + CELER_ASSERT(mucf_matid); + 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 + } + + //! \todo Correct track time update? Or should be done in Interactors? + 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) + { + case MucfMuonicMolecule::deuterium_deuterium: { + // Return DD interaction + DDMucfInteractor interact( + data, detail::DDChannelSelector()(rng), allocate_secondaries); + result = interact(rng); + break; + } + case MucfMuonicMolecule::deuterium_tritium: { + // Return DT interaction + DTMucfInteractor interact( + data, detail::DTChannelSelector()(rng), allocate_secondaries); + result = interact(rng); + break; + } + case MucfMuonicMolecule::tritium_tritium: { + // Return TT interaction + TTMucfInteractor interact( + data, detail::TTChannelSelector()(rng), allocate_secondaries); + result = interact(rng); + break; + } + default: + CELER_ASSERT_UNREACHABLE(); + } + + //! \todo Muon stripping: strip muon from muonic atom secondaries + // May be added as a separate discrete process in the stepping loop + + return result; +} + +//---------------------------------------------------------------------------// +} // namespace celeritas diff --git a/src/celeritas/mucf/executor/detail/DDChannelSelector.hh b/src/celeritas/mucf/executor/detail/DDChannelSelector.hh new file mode 100644 index 0000000000..05394c8019 --- /dev/null +++ b/src/celeritas/mucf/executor/detail/DDChannelSelector.hh @@ -0,0 +1,72 @@ +//------------------------------- -*- C++ -*- -------------------------------// +// Copyright Celeritas contributors: see top-level COPYRIGHT file for details +// SPDX-License-Identifier: (Apache-2.0 OR MIT) +//---------------------------------------------------------------------------// +//! \file celeritas/mucf/executor/detail/DDChannelSelector.hh +//---------------------------------------------------------------------------// +#pragma once + +#include "celeritas/mucf/interactor/DDMucfInteractor.hh" + +namespace celeritas +{ +namespace detail +{ +//---------------------------------------------------------------------------// +/*! + * Select final channel for muonic dd molecules. + * + * This selection already accounts for sticking, as that is one of the possible + * outcomes. + */ +class DDChannelSelector +{ + public: + //!@{ + //! \name Type aliases + using Channel = DDMucfInteractor::Channel; + //!@} + + public: + //! Construct with args; \todo Update documentation + inline CELER_FUNCTION DDChannelSelector(/* args */); + + // Select fusion channel to be used by the interactor + template + inline CELER_FUNCTION Channel operator()(Engine&); +}; + +//---------------------------------------------------------------------------// +// INLINE DEFINITIONS +//---------------------------------------------------------------------------// +/*! + * Construct with args. + * + * \todo Update documentation + */ +CELER_FUNCTION DDChannelSelector::DDChannelSelector(/* args */) +{ + CELER_NOT_IMPLEMENTED("Mucf dd fusion channel selection"); +} + +//---------------------------------------------------------------------------// +/*! + * Return a sampled channel to be used as input in the dd muCF interactor. + * + * \sa celeritas::DDMucfInteractor + */ +template +CELER_FUNCTION DDChannelSelector::Channel DDChannelSelector::operator()(Engine&) +{ + Channel result{Channel::size_}; + + //! \todo Implement + // Final channel selection already takes into account sticking. + + CELER_ENSURE(result < Channel::size_); + return result; +} + +//---------------------------------------------------------------------------// +} // namespace detail +} // namespace celeritas diff --git a/src/celeritas/mucf/executor/detail/DTChannelSelector.hh b/src/celeritas/mucf/executor/detail/DTChannelSelector.hh new file mode 100644 index 0000000000..47e097cdfd --- /dev/null +++ b/src/celeritas/mucf/executor/detail/DTChannelSelector.hh @@ -0,0 +1,72 @@ +//------------------------------- -*- C++ -*- -------------------------------// +// Copyright Celeritas contributors: see top-level COPYRIGHT file for details +// SPDX-License-Identifier: (Apache-2.0 OR MIT) +//---------------------------------------------------------------------------// +//! \file celeritas/mucf/executor/detail/DTChannelSelector.hh +//---------------------------------------------------------------------------// +#pragma once + +#include "celeritas/mucf/interactor/DTMucfInteractor.hh" + +namespace celeritas +{ +namespace detail +{ +//---------------------------------------------------------------------------// +/*! + * Select final channel for muonic dt molecules. + * + * This selection already accounts for sticking, as that is one of the possible + * outcomes. + */ +class DTChannelSelector +{ + public: + //!@{ + //! \name Type aliases + using Channel = DTMucfInteractor::Channel; + //!@} + + public: + //! Construct with args; \todo Update documentation + inline CELER_FUNCTION DTChannelSelector(/* args */); + + // Select fusion channel to be used by the interactor + template + inline CELER_FUNCTION Channel operator()(Engine&); +}; + +//---------------------------------------------------------------------------// +// 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. + * + * \sa celeritas::DTMucfInteractor + */ +template +CELER_FUNCTION DTChannelSelector::Channel DTChannelSelector::operator()(Engine&) +{ + Channel result{Channel::size_}; + + //! \todo Implement + // Final channel selection already takes into account sticking. + + CELER_ENSURE(result < Channel::size_); + return result; +} + +//---------------------------------------------------------------------------// +} // namespace detail +} // namespace celeritas diff --git a/src/celeritas/mucf/executor/detail/DTMixMuonicAtomSelector.hh b/src/celeritas/mucf/executor/detail/DTMixMuonicAtomSelector.hh new file mode 100644 index 0000000000..f193109e0c --- /dev/null +++ b/src/celeritas/mucf/executor/detail/DTMixMuonicAtomSelector.hh @@ -0,0 +1,61 @@ +//------------------------------- -*- C++ -*- -------------------------------// +// Copyright Celeritas contributors: see top-level COPYRIGHT file for details +// SPDX-License-Identifier: (Apache-2.0 OR MIT) +//---------------------------------------------------------------------------// +//! \file celeritas/mucf/executor/detail/DTMixMuonicAtomSelector.hh +//---------------------------------------------------------------------------// +#pragma once + +#include "corecel/cont/EnumArray.hh" +#include "celeritas/Types.hh" + +namespace celeritas +{ +namespace detail +{ +//---------------------------------------------------------------------------// +/*! + * Select a muonic atom given the material information. + */ +class DTMixMuonicAtomSelector +{ + public: + //! Construct with args; \todo Update documentation + inline CELER_FUNCTION DTMixMuonicAtomSelector(/* args */); + + // Select muonic atom + template + inline CELER_FUNCTION MucfMuonicAtom operator()(Engine&); +}; + +//---------------------------------------------------------------------------// +// INLINE DEFINITIONS +//---------------------------------------------------------------------------// +/*! + * Construct with args. + + * \todo Update documentation + */ +CELER_FUNCTION DTMixMuonicAtomSelector::DTMixMuonicAtomSelector(/* args */) +{ + CELER_NOT_IMPLEMENTED("Mucf muonic atom selection"); +} + +//---------------------------------------------------------------------------// +/*! + * Return selected muonic atom. + */ +template +CELER_FUNCTION MucfMuonicAtom DTMixMuonicAtomSelector::operator()(Engine&) +{ + MucfMuonicAtom result{MucfMuonicAtom::size_}; + + //! \todo Implement + + CELER_ENSURE(result < MucfMuonicAtom::size_); + return result; +} + +//---------------------------------------------------------------------------// +} // namespace detail +} // namespace celeritas diff --git a/src/celeritas/mucf/executor/detail/DTMixMuonicMoleculeSelector.hh b/src/celeritas/mucf/executor/detail/DTMixMuonicMoleculeSelector.hh new file mode 100644 index 0000000000..cc4c7bd30b --- /dev/null +++ b/src/celeritas/mucf/executor/detail/DTMixMuonicMoleculeSelector.hh @@ -0,0 +1,67 @@ +//------------------------------- -*- 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/MuonicAtomSpinSelector.hh b/src/celeritas/mucf/executor/detail/MuonicAtomSpinSelector.hh new file mode 100644 index 0000000000..6d22817345 --- /dev/null +++ b/src/celeritas/mucf/executor/detail/MuonicAtomSpinSelector.hh @@ -0,0 +1,64 @@ +//------------------------------- -*- C++ -*- -------------------------------// +// Copyright Celeritas contributors: see top-level COPYRIGHT file for details +// SPDX-License-Identifier: (Apache-2.0 OR MIT) +//---------------------------------------------------------------------------// +//! \file celeritas/mucf/executor/detail/MuonicAtomSpinSelector.hh +//---------------------------------------------------------------------------// +#pragma once + +#include "corecel/Assert.hh" +#include "celeritas/Types.hh" + +namespace celeritas +{ +namespace detail +{ +//---------------------------------------------------------------------------// +/*! + * Select muonic atom spin, in units of \f$ \frac{\hbar}{2} \f$. + * + * Applicable to \f$ (p)_\mu \f$, \f$ (d)_\mu \f$, and \f$ (t)_\mu \f$. + */ +class MuonicAtomSpinSelector +{ + public: + // Construct with muonic atom type + inline CELER_FUNCTION MuonicAtomSpinSelector(MucfMuonicAtom atom); + + // Sample and return a spin value in units of hbar / 2 + template + inline CELER_FUNCTION size_type operator()(Engine&); + + private: + MucfMuonicAtom atom_; + //! \todo Add constant atom spin sampling rejection fractions +}; + +//---------------------------------------------------------------------------// +// INLINE DEFINITIONS +//---------------------------------------------------------------------------// +/*! + * Construct with muonic atom. + */ +CELER_FUNCTION +MuonicAtomSpinSelector::MuonicAtomSpinSelector(MucfMuonicAtom atom) + : atom_(atom) +{ + CELER_EXPECT(atom_ < MucfMuonicAtom::size_); + CELER_NOT_IMPLEMENTED("Mucf muonic atom spin selection"); +} + +//---------------------------------------------------------------------------// +/*! + * Return selected spin, in units of \f$ \hbar / 2 \f$. + */ +template +CELER_FUNCTION size_type MuonicAtomSpinSelector::operator()(Engine&) +{ + //! \todo switch on atom_ + CELER_ASSERT_UNREACHABLE(); +} + +//---------------------------------------------------------------------------// +} // namespace detail +} // namespace celeritas diff --git a/src/celeritas/mucf/executor/detail/MuonicMoleculeSpinSelector.hh b/src/celeritas/mucf/executor/detail/MuonicMoleculeSpinSelector.hh new file mode 100644 index 0000000000..67ae8d4e29 --- /dev/null +++ b/src/celeritas/mucf/executor/detail/MuonicMoleculeSpinSelector.hh @@ -0,0 +1,64 @@ +//------------------------------- -*- 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 new file mode 100644 index 0000000000..1f7c0b40ea --- /dev/null +++ b/src/celeritas/mucf/executor/detail/TTChannelSelector.hh @@ -0,0 +1,73 @@ +//------------------------------- -*- C++ -*- -------------------------------// +// Copyright Celeritas contributors: see top-level COPYRIGHT file for details +// SPDX-License-Identifier: (Apache-2.0 OR MIT) +//---------------------------------------------------------------------------// +//! \file celeritas/mucf/executor/detail/TTChannelSelectionHelper.hh +//---------------------------------------------------------------------------// +#pragma once + +#include "celeritas/mucf/interactor/TTMucfInteractor.hh" + +namespace celeritas +{ +namespace detail +{ +//---------------------------------------------------------------------------// +/*! + * Select final channel for muonic tt molecules. + * + * This selection already accounts for sticking, as that is one of the possible + * outcomes. + */ +class TTChannelSelector +{ + public: + //!@{ + //! \name Type aliases + using Channel = TTMucfInteractor::Channel; + //!@} + + public: + //! Construct with args; \todo Update documentation + inline CELER_FUNCTION TTChannelSelector(/* args */); + + // Select fusion channel to be used by the interactor + template + inline CELER_FUNCTION Channel operator()(Engine&); +}; + +//---------------------------------------------------------------------------// +// 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. + * + * \sa celeritas::TTMucfInteractor + */ +template +CELER_FUNCTION TTChannelSelector::Channel TTChannelSelector::operator()(Engine&) +{ + Channel result{Channel::size_}; + + //! \todo Implement + // Final channel selection already takes into account sticking. + + CELER_ENSURE(result < Channel::size_); + return result; +} + +//---------------------------------------------------------------------------// +} // namespace detail +} // namespace celeritas diff --git a/src/celeritas/mucf/interactor/DDMucfInteractor.hh b/src/celeritas/mucf/interactor/DDMucfInteractor.hh new file mode 100644 index 0000000000..b75fa71c02 --- /dev/null +++ b/src/celeritas/mucf/interactor/DDMucfInteractor.hh @@ -0,0 +1,128 @@ +//------------------------------- -*- C++ -*- -------------------------------// +// Copyright Celeritas contributors: see top-level COPYRIGHT file for details +// SPDX-License-Identifier: (Apache-2.0 OR MIT) +//---------------------------------------------------------------------------// +//! \file celeritas/mucf/interactor/DDMucfInteractor.hh +//---------------------------------------------------------------------------// +#pragma once + +#include "corecel/Macros.hh" +#include "corecel/data/StackAllocator.hh" +#include "celeritas/mat/ElementView.hh" +#include "celeritas/mat/MaterialView.hh" +#include "celeritas/mucf/data/DTMixMucfData.hh" +#include "celeritas/phys/Interaction.hh" +#include "celeritas/phys/ParticleTrackView.hh" +#include "celeritas/phys/Secondary.hh" + +namespace celeritas +{ +//---------------------------------------------------------------------------// +/*! + * Muon-catalyzed fusion of \f$ (dd)_\mu \f$ molecules. + * + * Fusion channels: + * - \f$ ^3\text{He} + \mu + n \f$ + * - \f$ (^3\text{He})_\mu + n \f$ + * - \f$ ^3\text{H} + \mu + p \f$ + */ +class DDMucfInteractor +{ + public: + //! \todo Implement muonichydrogen3_proton (\f$ (^3\text{H})_\mu + p \f$) + enum class Channel + { + helium3_muon_neutron, //!< \f$ ^3\text{He} + \mu + n \f$ + muonichelium3_neutron, //!< \f$ (^3\text{He})_\mu + n \f$ + hydrogen3_muon_proton, //!< \f$ ^3\text{H} + \mu + p \f$ + size_ + }; + + // Construct from shared and state data + inline CELER_FUNCTION + DDMucfInteractor(NativeCRef const& data, + Channel channel, + StackAllocator& allocate); + + // Sample an interaction with the given RNG + template + inline CELER_FUNCTION Interaction operator()(Engine& rng); + + private: + // Shared constant physics properties + NativeCRef const& data_; + // Selected fusion channel + Channel channel_{Channel::size_}; + // Allocate space for secondary particles + StackAllocator& allocate_; + // Number of secondaries per channel + EnumArray num_secondaries_{ + 3, // helium3_muon_neutron + 2, // muonichelium3_neutron + 3 // hydrogen3_muon_proton + }; + + // Sample Interaction secondaries + template + inline CELER_FUNCTION Span + sample_secondaries(Secondary* secondaries /*, other args */, Engine&); +}; + +//---------------------------------------------------------------------------// +// INLINE DEFINITIONS +//---------------------------------------------------------------------------// +/*! + * Construct with shared and state data. + */ +CELER_FUNCTION +DDMucfInteractor::DDMucfInteractor(NativeCRef const& data, + Channel const channel, + StackAllocator& allocate) + : data_(data), channel_(channel), allocate_(allocate) +{ + CELER_EXPECT(data_); + CELER_EXPECT(channel < Channel::size_); +} + +//---------------------------------------------------------------------------// +/*! + * Sample a dd muonic molecule fusion. + */ +template +CELER_FUNCTION Interaction DDMucfInteractor::operator()(Engine& rng) +{ + // Allocate space for the final fusion channel + Secondary* secondaries = allocate_(num_secondaries_[channel_]); + if (secondaries == nullptr) + { + // Failed to allocate space for secondaries + return Interaction::from_failure(); + } + + // Kill primary and generate secondaries + Interaction result = Interaction::from_absorption(); + result.secondaries = this->sample_secondaries(secondaries, rng); + + return result; +} + +//---------------------------------------------------------------------------// +/*! + * Sample the secondaries of the selected channel. + * + * Since secondaries come from an at rest interaction, their final state is + * a simple combination of random direction + momentum conservation + */ +template +CELER_FUNCTION Span +DDMucfInteractor::sample_secondaries(Secondary* secondaries /*, other args */, + Engine&) +{ + // TODO: switch on channel_ + CELER_ASSERT_UNREACHABLE(); + + return Span{secondaries, num_secondaries_[channel_]}; +} + +//---------------------------------------------------------------------------// +} // namespace celeritas diff --git a/src/celeritas/mucf/interactor/DTMucfInteractor.hh b/src/celeritas/mucf/interactor/DTMucfInteractor.hh new file mode 100644 index 0000000000..568bf3ba0e --- /dev/null +++ b/src/celeritas/mucf/interactor/DTMucfInteractor.hh @@ -0,0 +1,123 @@ +//------------------------------- -*- C++ -*- -------------------------------// +// Copyright Celeritas contributors: see top-level COPYRIGHT file for details +// SPDX-License-Identifier: (Apache-2.0 OR MIT) +//---------------------------------------------------------------------------// +//! \file celeritas/mucf/interactor/DTMucfInteractor.hh +//---------------------------------------------------------------------------// +#pragma once + +#include "corecel/Macros.hh" +#include "corecel/data/StackAllocator.hh" +#include "celeritas/mat/ElementView.hh" +#include "celeritas/mat/MaterialView.hh" +#include "celeritas/mucf/data/DTMixMucfData.hh" +#include "celeritas/phys/Interaction.hh" +#include "celeritas/phys/ParticleTrackView.hh" +#include "celeritas/phys/Secondary.hh" + +namespace celeritas +{ +//---------------------------------------------------------------------------// +/*! + * Muon-catalyzed fusion of \f$ (dt)_\mu \f$ molecules. + * + * Fusion channels: + * - \f$ \alpha + \mu + n \f$ + * - \f$ (\alpha)_\mu + n \f$ + */ +class DTMucfInteractor +{ + public: + enum class Channel + { + alpha_muon_neutron, //!< \f$ \alpha + \mu + n \f$ + muonicalpha_neutron, //!< \f$ (\alpha)_\mu + n \f$ + size_ + }; + + // Construct from shared and state data + inline CELER_FUNCTION + DTMucfInteractor(NativeCRef const& data, + Channel channel, + StackAllocator& allocate); + + // Sample an interaction with the given RNG + template + inline CELER_FUNCTION Interaction operator()(Engine& rng); + + private: + // Shared constant physics properties + NativeCRef const& data_; + // Selected fusion channel + Channel channel_{Channel::size_}; + // Allocate space for secondary particles + StackAllocator& allocate_; + // Number of secondaries per channel + EnumArray num_secondaries_{ + 3, // alpha_muon_neutron + 2 // muonicalpha_neutron + }; + + // Sample Interaction secondaries + template + inline CELER_FUNCTION Span + sample_secondaries(Secondary* secondaries /*, other args */, Engine&); +}; + +//---------------------------------------------------------------------------// +// INLINE DEFINITIONS +//---------------------------------------------------------------------------// +/*! + * Construct with shared and state data. + */ +CELER_FUNCTION +DTMucfInteractor::DTMucfInteractor(NativeCRef const& data, + Channel const channel, + StackAllocator& allocate) + : data_(data), channel_(channel), allocate_(allocate) +{ + CELER_EXPECT(data_); + CELER_EXPECT(channel_ < Channel::size_); +} + +//---------------------------------------------------------------------------// +/*! + * Sample a dt muonic molecule fusion. + */ +template +CELER_FUNCTION Interaction DTMucfInteractor::operator()(Engine& rng) +{ + // Allocate space for the final fusion channel + Secondary* secondaries = allocate_(num_secondaries_[channel_]); + if (secondaries == nullptr) + { + // Failed to allocate space for secondaries + return Interaction::from_failure(); + } + + // Kill primary and generate secondaries + Interaction result = Interaction::from_absorption(); + result.secondaries = this->sample_secondaries(secondaries, rng); + + return result; +} + +//---------------------------------------------------------------------------// +/*! + * Sample the secondaries of the selected channel. + * + * Since secondaries come from an at rest interaction, their final state is + * a simple combination of random direction + momentum conservation + */ +template +CELER_FUNCTION Span +DTMucfInteractor::sample_secondaries(Secondary* secondaries /*, other args */, + Engine&) +{ + // TODO: switch on channel_ + CELER_ASSERT_UNREACHABLE(); + return Span{secondaries, num_secondaries_[channel_]}; +} + +//---------------------------------------------------------------------------// +} // namespace celeritas diff --git a/src/celeritas/mucf/interactor/TTMucfInteractor.hh b/src/celeritas/mucf/interactor/TTMucfInteractor.hh new file mode 100644 index 0000000000..d0239fc561 --- /dev/null +++ b/src/celeritas/mucf/interactor/TTMucfInteractor.hh @@ -0,0 +1,124 @@ +//------------------------------- -*- C++ -*- -------------------------------// +// Copyright Celeritas contributors: see top-level COPYRIGHT file for details +// SPDX-License-Identifier: (Apache-2.0 OR MIT) +//---------------------------------------------------------------------------// +//! \file celeritas/mucf/interactor/TTMucfInteractor.hh +//---------------------------------------------------------------------------// +#pragma once + +#include "corecel/Macros.hh" +#include "corecel/data/StackAllocator.hh" +#include "celeritas/mat/ElementView.hh" +#include "celeritas/mat/MaterialView.hh" +#include "celeritas/mucf/data/DTMixMucfData.hh" +#include "celeritas/phys/Interaction.hh" +#include "celeritas/phys/ParticleTrackView.hh" +#include "celeritas/phys/Secondary.hh" + +namespace celeritas +{ +//---------------------------------------------------------------------------// +/*! + * Muon-catalyzed fusion of \f$ (tt)_\mu \f$ molecules. + * + * Fusion channels: + * - \f$ \alpha + \mu + n + n \f$ + * - \f$ (\alpha)_\mu + n + n \f$ + */ +class TTMucfInteractor +{ + public: + enum class Channel + { + alpha_muon_neutron_neutron, //!< \f$ \alpha + \mu + n + n \f$ + muonicalpha_neutron_neutron, //!< \f$ (\alpha)_\mu + n + n \f$ + size_ + }; + + // Construct from shared and state data + inline CELER_FUNCTION + TTMucfInteractor(NativeCRef const& data, + Channel channel, + StackAllocator& allocate); + + // Sample an interaction with the given RNG + template + inline CELER_FUNCTION Interaction operator()(Engine& rng); + + private: + // Shared constant physics properties + NativeCRef const& data_; + // Selected fusion channel + Channel channel_{Channel::size_}; + // Allocate space for secondary particles + StackAllocator& allocate_; + // Number of secondaries per channel + EnumArray num_secondaries_{ + 4, // alpha_muon_neutron_neutron + 3 // muonicalpha_neutron_neutron + }; + + // Sample Interaction secondaries + template + inline CELER_FUNCTION Span + sample_secondaries(Secondary* secondaries /*, other args */, Engine&); +}; + +//---------------------------------------------------------------------------// +// INLINE DEFINITIONS +//---------------------------------------------------------------------------// +/*! + * Construct with shared and state data. + */ +CELER_FUNCTION +TTMucfInteractor::TTMucfInteractor(NativeCRef const& data, + Channel const channel, + StackAllocator& allocate) + : data_(data), channel_(channel), allocate_(allocate) +{ + CELER_EXPECT(data_); + CELER_EXPECT(channel_ < Channel::size_); +} + +//---------------------------------------------------------------------------// +/*! + * Sample a tt muonic molecule fusion. + */ +template +CELER_FUNCTION Interaction TTMucfInteractor::operator()(Engine& rng) +{ + // Allocate space for the final fusion channel + Secondary* secondaries = allocate_(num_secondaries_[channel_]); + if (secondaries == nullptr) + { + // Failed to allocate space for secondaries + return Interaction::from_failure(); + } + + // Kill primary and generate secondaries + Interaction result = Interaction::from_absorption(); + result.secondaries = this->sample_secondaries(secondaries, rng); + + return result; +} + +//---------------------------------------------------------------------------// +/*! + * Sample the secondaries of the selected channel. + * + * Since secondaries come from an at rest interaction, their final state is + * a simple combination of random direction + momentum conservation + */ +template +CELER_FUNCTION Span +TTMucfInteractor::sample_secondaries(Secondary* secondaries /*, other args */, + Engine&) +{ + // TODO: switch on channel_ + CELER_ASSERT_UNREACHABLE(); + + return Span{secondaries, num_secondaries_[channel_]}; +} + +//---------------------------------------------------------------------------// +} // namespace celeritas diff --git a/src/celeritas/mucf/model/DTMixMucfModel.cc b/src/celeritas/mucf/model/DTMixMucfModel.cc new file mode 100644 index 0000000000..d4cee2333a --- /dev/null +++ b/src/celeritas/mucf/model/DTMixMucfModel.cc @@ -0,0 +1,179 @@ +//------------------------------- -*- C++ -*- -------------------------------// +// Copyright Celeritas contributors: see top-level COPYRIGHT file for details +// SPDX-License-Identifier: (Apache-2.0 OR MIT) +//---------------------------------------------------------------------------// +//! \file celeritas/mucf/model/DTMixMucfModel.cc +//---------------------------------------------------------------------------// +#include "DTMixMucfModel.hh" + +#include + +#include "corecel/OpaqueIdIO.hh" +#include "corecel/inp/Grid.hh" +#include "celeritas/global/ActionLauncher.hh" +#include "celeritas/global/TrackExecutor.hh" +#include "celeritas/grid/NonuniformGridBuilder.hh" +#include "celeritas/inp/MucfPhysics.hh" +#include "celeritas/mat/MaterialParams.hh" +#include "celeritas/mucf/executor/DTMixMucfExecutor.hh" // IWYU pragma: associated +#include "celeritas/phys/InteractionApplier.hh" // IWYU pragma: associated +#include "celeritas/phys/PDGNumber.hh" +#include "celeritas/phys/ParticleParams.hh" + +#include "detail/MucfMaterialInserter.hh" + +namespace celeritas +{ +namespace +{ +//---------------------------------------------------------------------------// +/*! + * Assign particle IDs from \c ParticleParams . + */ +static MucfParticleIds from_params(ParticleParams const& particles) +{ + using PairStrPdg = std::pair; + std::vector missing; + MucfParticleIds result; + +#define MP_ADD(MEMBER) \ + result.MEMBER = particles.find(pdg::MEMBER()); \ + if (!result.MEMBER) \ + { \ + missing.push_back({#MEMBER, pdg::MEMBER()}); \ + } + + MP_ADD(mu_minus); + MP_ADD(neutron); + MP_ADD(proton); + MP_ADD(alpha); + MP_ADD(he3); + MP_ADD(muonic_deuteron); + MP_ADD(muonic_triton); + + //! \todo Decide whether to implement these PDGs in PDGNumber.hh +#if 0 + MP_ADD(muonic_hydrogen); + MP_ADD(muonic_alpha); + MP_ADD(muonic_he3); +#endif + + CELER_VALIDATE(missing.empty(), + << "missing particles required for muon-catalyzed fusion: " + << join_stream(missing.begin(), + missing.end(), + ", ", + [](std::ostream& os, PairStrPdg const& p) { + os << p.first << " (PDG " + << p.second.unchecked_get() << ')'; + })); + return result; + +#undef MP_ADD +} +//---------------------------------------------------------------------------// +} // namespace + +//---------------------------------------------------------------------------// +/*! + * Construct from model ID and other necessary data. + * + * Most of the muon-catalyzed fusion data is static throughout the simulation, + * as it is only material-dependent (DT mixture and temperature). Therefore, + * most grids can be host-only and used to calculate final values, which are + * then cached and copied to device. The exception to this is the muon energy + * CDF grid, needed to sample the final state of the outgoing muon after a muCF + * interaction. + * + * \todo Correctly update \c ImportProcessClass and \c ImportModelClass . These + * operate under the assumption that there is a one-to-one equivalente between + * Geant4 and Celeritas. But for muCF, everything is done via one + * process/model/executor in Celeritas, whereas in Geant4 atom formation, spin + * flip, atom transfer, etc., are are all separate processes. + */ +DTMixMucfModel::DTMixMucfModel(ActionId id, + ParticleParams const& particles, + MaterialParams const& materials) + : StaticConcreteAction( + id, + "dt-mucf", + R"(interact by muon forming and fusing a dd, dt, or tt muonic molecule)") +{ + CELER_EXPECT(id); + + // Initialize muCF physics input data + inp::MucfPhysics inp_data = inp::MucfPhysics::from_default(); + CELER_EXPECT(inp_data); + + HostVal host_data; + host_data.particles = from_params(particles); + + // Copy muon energy CDF data using NonuniformGridBuilder + NonuniformGridBuilder build_grid_record{&host_data.reals}; + host_data.muon_energy_cdf = build_grid_record(inp_data.muon_energy_cdf); + + // Calculate and cache quantities for all materials with dt mixtures + detail::MucfMaterialInserter insert(&host_data); + for (auto const& matid : range(materials.num_materials())) + { + auto const& mat_view = materials.get(PhysMatId{matid}); + if (insert(mat_view)) + { + CELER_LOG(debug) << "Added material ID " << mat_view.material_id() + << " as a muCF d-t mixture"; + } + } + + // Copy to device + data_ = CollectionMirror{std::move(host_data)}; + CELER_ENSURE(this->data_); +} + +//---------------------------------------------------------------------------// +/*! + * Particle types and energy ranges that this model applies to. + */ +auto DTMixMucfModel::applicability() const -> SetApplicability +{ + Applicability applic; + applic.particle = this->host_ref().particles.mu_minus; + // At-rest model + applic.lower = zero_quantity(); + applic.upper = zero_quantity(); + + return {applic}; +} + +//---------------------------------------------------------------------------// +/*! + * At-rest model does not require microscopic cross sections. + */ +auto DTMixMucfModel::micro_xs(Applicability) const -> XsTable +{ + return {}; +} + +//---------------------------------------------------------------------------// +/*! + * Interact with host data. + */ +void DTMixMucfModel::step(CoreParams const& params, CoreStateHost& state) const +{ + auto execute = make_action_track_executor( + params.ptr(), + state.ptr(), + this->action_id(), + InteractionApplier{DTMixMucfExecutor{this->host_ref()}}); + return launch_action(*this, params, state, execute); +} + +//---------------------------------------------------------------------------// +#if !CELER_USE_DEVICE +void DTMixMucfModel::step(CoreParams const&, CoreStateDevice&) const +{ + CELER_NOT_CONFIGURED("CUDA OR HIP"); +} +#endif + +//---------------------------------------------------------------------------// +} // namespace celeritas diff --git a/src/celeritas/mucf/model/DTMixMucfModel.cu b/src/celeritas/mucf/model/DTMixMucfModel.cu new file mode 100644 index 0000000000..4ac694e935 --- /dev/null +++ b/src/celeritas/mucf/model/DTMixMucfModel.cu @@ -0,0 +1,33 @@ +//------------------------------- -*- C++ -*- -------------------------------// +// Copyright Celeritas contributors: see top-level COPYRIGHT file for details +// SPDX-License-Identifier: (Apache-2.0 OR MIT) +//---------------------------------------------------------------------------// +//! \file celeritas/mucf/model/DTMixMucfModel.cu +//---------------------------------------------------------------------------// +#include "DTMixMucfModel.hh" + +#include "celeritas/global/ActionLauncher.device.hh" +#include "celeritas/global/CoreParams.hh" +#include "celeritas/global/CoreState.hh" +#include "celeritas/global/TrackExecutor.hh" +#include "celeritas/mucf/executor/DTMixMucfExecutor.hh" +#include "celeritas/phys/InteractionApplier.hh" + +namespace celeritas +{ +//---------------------------------------------------------------------------// +/*! + * Interact with device data. + */ +void DTMixMucfModel::step(CoreParams const& params, CoreStateDevice& state) const +{ + auto execute = make_action_track_executor( + params.ptr(), + state.ptr(), + this->action_id(), + InteractionApplier{DTMixMucfExecutor{this->device_ref()}}); + static ActionLauncher const launch_kernel(*this); + launch_kernel(*this, params, state, execute); +} +//---------------------------------------------------------------------------// +} // namespace celeritas diff --git a/src/celeritas/mucf/model/DTMixMucfModel.hh b/src/celeritas/mucf/model/DTMixMucfModel.hh new file mode 100644 index 0000000000..893d6c2e5e --- /dev/null +++ b/src/celeritas/mucf/model/DTMixMucfModel.hh @@ -0,0 +1,74 @@ +//------------------------------- -*- C++ -*- -------------------------------// +// Copyright Celeritas contributors: see top-level COPYRIGHT file for details +// SPDX-License-Identifier: (Apache-2.0 OR MIT) +//---------------------------------------------------------------------------// +//! \file celeritas/mucf/model/DTMixMucfModel.hh +//---------------------------------------------------------------------------// +#pragma once + +#include "celeritas/mucf/data/DTMixMucfData.hh" +#include "celeritas/phys/Model.hh" + +namespace celeritas +{ +class MaterialParams; +class ParticleParams; + +//---------------------------------------------------------------------------// +/*! + * Muon-catalyzed fusion model for dd, dt, and tt molecules. + * + * In this model the executor performs the full muon-catalyzed fusion workflow. + * It forms a muonic d or t atom, samples which muonic molecule will be + * produced, selects the channel, and calls the appropriate interactor. + * + * The full set of "actions" is as follows, and in this ordering: + * - Define muon decay time to compete with the rest of the execution + * - Form muonic atom and select its spin + * - May execute atom spin flip or atom transfer + * - Form muonic molecule and select its spin + * - Calculate mean cycle time (time it takes from atom formation to fusion) + * - Confirm if fusion happens or the if the muon should decay + * - Call appropriate Interactor: Muon decay, or one of the muCF interactors + * + * \note This is an at-rest model. + */ +class DTMixMucfModel final : public Model, public StaticConcreteAction +{ + public: + //!@{ + //! \name Type aliases + using HostRef = HostCRef; + using DeviceRef = DeviceCRef; + //!@} + + // Construct from model ID and other necessary data + DTMixMucfModel(ActionId id, + ParticleParams const& particles, + MaterialParams const& materials); + + // Particle types and energy ranges that this model applies to + SetApplicability applicability() const final; + + // Get the microscopic cross sections for the given particle and material + XsTable micro_xs(Applicability) const final; + + // Apply the interaction kernel on host + void step(CoreParams const&, CoreStateHost&) const final; + + // Apply the interaction kernel on device + void step(CoreParams const&, CoreStateDevice&) const final; + + //!@{ + //! Access model data + HostRef const& host_ref() const { return data_.host_ref(); } + DeviceRef const& device_ref() const { return data_.device_ref(); } + //!@} + + private: + // Host/device storage and reference + CollectionMirror data_; +}; + +//---------------------------------------------------------------------------// +} // namespace celeritas diff --git a/src/celeritas/mucf/model/detail/MucfMaterialInserter.cc b/src/celeritas/mucf/model/detail/MucfMaterialInserter.cc new file mode 100644 index 0000000000..f859297ab9 --- /dev/null +++ b/src/celeritas/mucf/model/detail/MucfMaterialInserter.cc @@ -0,0 +1,243 @@ +//------------------------------- -*- C++ -*- -------------------------------// +// Copyright Celeritas contributors: see top-level COPYRIGHT file for details +// SPDX-License-Identifier: (Apache-2.0 OR MIT) +//---------------------------------------------------------------------------// +//! \file celeritas/mucf/model/detail/MucfMaterialInserter.cc +//---------------------------------------------------------------------------// +#include "MucfMaterialInserter.hh" + +#include "corecel/Assert.hh" + +namespace celeritas +{ +namespace detail +{ +//---------------------------------------------------------------------------// +/*! + * Construct with \c DTMixMucfModel model data. + */ +MucfMaterialInserter::MucfMaterialInserter(HostVal* host_data) + : mucfmatid_to_matid_(&host_data->mucfmatid_to_matid) + , cycle_times_(&host_data->cycle_times) +{ + CELER_EXPECT(host_data); +} + +//---------------------------------------------------------------------------// +/*! + * Insert material information if applicable. + * + * Calculates and caches material-dependent properties needed by the + * \c DTMixMucfModel . If the material does not contain deuterium and/or + * tritium the operator will return false. + */ +bool MucfMaterialInserter::operator()(MaterialView const& material) +{ + for (auto elcompid : range(material.num_elements())) + { + auto const& element_view + = material.element_record(ElementComponentId{elcompid}); + if (element_view.atomic_number() != AtomicNumber{1}) + { + // Skip non-hydrogen elements + continue; + } + + // Found hydrogen; Check for d and t isotopes + IsotopeChecker has_isotope{false, false}; + for (auto el_comp : range(element_view.num_isotopes())) + { + auto iso_view + = element_view.isotope_record(IsotopeComponentId{el_comp}); + auto mass = iso_view.atomic_mass_number(); + if (mass == AtomicMassNumber{1}) + { + // Skip protium + continue; + } + + auto const atom = this->from_mass_number(mass); + if (atom < MucfMuonicAtom::size_) + { + // Mark d or t isotope as found + has_isotope[atom] = true; + } + } + + if (!has_isotope[MucfMuonicAtom::deuterium] + && !has_isotope[MucfMuonicAtom::tritium]) + { + // No deuterium or tritium found; skip material + return false; + } + + // Temporary data needed to calculate model data, such as cycle times + lhd_densities_ = this->calc_lhd_densities(element_view); + equilibrium_densities_ = this->calc_equilibrium_densities(element_view); + + // Calculate and insert muCF material data into model data + mucfmatid_to_matid_.push_back(material.material_id()); + cycle_times_.push_back( + this->calc_cycle_times(element_view, has_isotope)); + //! \todo Store mean atom spin flip and transfer times + } + return true; +} + +//---------------------------------------------------------------------------// +/*! + * Return \c MucfMuonicAtom from a given atomic mass number. + */ +MucfMuonicAtom MucfMaterialInserter::from_mass_number(AtomicMassNumber mass) +{ + if (mass == AtomicMassNumber{2}) + { + return MucfMuonicAtom::deuterium; + } + if (mass == AtomicMassNumber{3}) + { + return MucfMuonicAtom::tritium; + } + return MucfMuonicAtom::size_; +} + +//---------------------------------------------------------------------------// +/*! + * Convert dt mixture densities to units of liquid hydrogen density. + * + * Used during cycle time calculations. + */ +MucfMaterialInserter::LhdArray +MucfMaterialInserter::calc_lhd_densities(ElementView const&) +{ + LhdArray result; + + //! \todo Implement + + return result; +} + +//---------------------------------------------------------------------------// +/*! + * Calculate dt mixture densities after reaching thermodynamical + * equilibrium. + * + * Used during cycle time calculations. + */ +MucfMaterialInserter::EquilibriumArray +MucfMaterialInserter::calc_equilibrium_densities(ElementView const&) +{ + EquilibriumArray result; + + //! \todo Implement + + return result; +} + +//---------------------------------------------------------------------------// +/*! + * Calculate fusion mean cycle times. + * + * This is designed to work with the user's material definition being either: + * - Single element, multiple isotopes (H element, with H, d, and t isotopes); + * or + * - Multiple elements, single isotope each (separate H, d, and t elements). + */ +MucfMaterialInserter::CycleTimesArray +MucfMaterialInserter::calc_cycle_times(ElementView const& element, + IsotopeChecker const& has_isotope) +{ + CycleTimesArray result; + for (auto el_comp : range(element.num_isotopes())) + { + auto iso_view = element.isotope_record(IsotopeComponentId{el_comp}); + + // Select possible muonic atom based on the isotope/element mass number + auto atom = this->from_mass_number(iso_view.atomic_mass_number()); + switch (atom) + { + // Calculate cycle times for dd molecules + case MucfMuonicAtom::deuterium: { + result[MucfMuonicMolecule::deuterium_deuterium] + = this->calc_dd_cycle(element); + if (has_isotope[MucfMuonicAtom::tritium]) + { + // Calculate cycle times for dt molecules + result[MucfMuonicMolecule::deuterium_tritium] + = this->calc_dt_cycle(element); + } + break; + } + // Calculate cycle times for tt molecules + case MucfMuonicAtom::tritium: { + result[MucfMuonicMolecule::tritium_tritium] + = this->calc_tt_cycle(element); + break; + } + default: + CELER_ASSERT_UNREACHABLE(); + } + } + return result; +} + +//---------------------------------------------------------------------------// +/*! + * Calculate dd muonic molecules cycle times from material properties and grid + * data. + * + * Cycle times for dd molecules come from F = 0 and F = 1 spin states. + */ +MucfMaterialInserter::MoleculeCycles +MucfMaterialInserter::calc_dd_cycle(ElementView const&) +{ + MoleculeCycles result; + + //! \todo Implement + + // Reactive states are F = 0 and F = 1 + CELER_ENSURE(result[0] >= 0 && result[1] >= 0); + return result; +} + +//---------------------------------------------------------------------------// +/*! + * Calculate dt muonic molecules cycle times from material properties and grid + * data. + * + * Cycle times for dt molecules come from F = 1/2 and F = 3/2 spin states. + */ +MucfMaterialInserter::MoleculeCycles +MucfMaterialInserter::calc_dt_cycle(ElementView const&) +{ + MoleculeCycles result; + + //! \todo Implement + + // Reactive states are F = 1/2 and F = 3/2 + CELER_ENSURE(result[0] >= 0 && result[1] >= 0); + return result; +} + +//---------------------------------------------------------------------------// +/*! + * Calculate tt muonic molecules cycle times from material properties and grid + * data. + * + * Cycle times for tt molecules come only from the F = 1/2 spin state. + */ +MucfMaterialInserter::MoleculeCycles +MucfMaterialInserter::calc_tt_cycle(ElementView const&) +{ + MoleculeCycles result; + + //! \todo Implement + + // Only F = 1/2 is reactive + CELER_ENSURE(result[0] >= 0 && result[1] == 0); + return result; +} + +//---------------------------------------------------------------------------// +} // namespace detail +} // namespace celeritas diff --git a/src/celeritas/mucf/model/detail/MucfMaterialInserter.hh b/src/celeritas/mucf/model/detail/MucfMaterialInserter.hh new file mode 100644 index 0000000000..ad7c307b9b --- /dev/null +++ b/src/celeritas/mucf/model/detail/MucfMaterialInserter.hh @@ -0,0 +1,77 @@ +//------------------------------- -*- C++ -*- -------------------------------// +// Copyright Celeritas contributors: see top-level COPYRIGHT file for details +// SPDX-License-Identifier: (Apache-2.0 OR MIT) +//---------------------------------------------------------------------------// +//! \file celeritas/mucf/model/detail/MucfMaterialInserter.hh +//---------------------------------------------------------------------------// +#pragma once + +#include "corecel/data/CollectionBuilder.hh" +#include "celeritas/mat/MaterialView.hh" +#include "celeritas/mucf/Types.hh" +#include "celeritas/mucf/data/DTMixMucfData.hh" + +namespace celeritas +{ +namespace detail +{ +//---------------------------------------------------------------------------// +/*! + * Helper class to calculate and insert muCF material-dependent data into + * \c DTMixMucfData . + */ +class MucfMaterialInserter +{ + public: + // Construct with muCF host data + explicit MucfMaterialInserter(HostVal* host_data); + + // Insert material if it is a valid d-t mixture + bool operator()(MaterialView const& material); + + private: + //// DATA //// + + using MoleculeCycles = Array; + using CycleTimesArray = EnumArray; + using LhdArray = EnumArray; + using EquilibriumArray = EnumArray; + using AtomicMassNumber = AtomicNumber; + using MucfIsotope = MucfMuonicAtom; + using IsotopeChecker = EnumArray; + + // DTMixMucfModel host data references populated by operator() + CollectionBuilder mucfmatid_to_matid_; + CollectionBuilder cycle_times_; + // Temporary quantities needed for calculating the model data + LhdArray lhd_densities_; + EquilibriumArray equilibrium_densities_; + + //// HELPER FUNCTIONS //// + + // Return muonic atom from given atomic mass number + MucfMuonicAtom from_mass_number(AtomicMassNumber mass); + + // Calculate dt mixture densities in units of liquid hydrogen density + LhdArray calc_lhd_densities(ElementView const&); + + // Calculate thermal equilibrium densities + EquilibriumArray calc_equilibrium_densities(ElementView const&); + + // Calculate mean fusion cycle times for all reactive muonic molecules + CycleTimesArray calc_cycle_times(ElementView const& element, + IsotopeChecker const& has_isotope); + + // Calculate mean fusion cycle times for dd muonic molecules + Array calc_dd_cycle(ElementView const&); + + // Calculate mean fusion cycle times for dt muonic molecules + Array calc_dt_cycle(ElementView const&); + + // Calculate mean fusion cycle times for tt muonic molecules + Array calc_tt_cycle(ElementView const&); +}; + +//---------------------------------------------------------------------------// +} // namespace detail +} // namespace celeritas diff --git a/src/celeritas/mucf/process/MucfProcess.cc b/src/celeritas/mucf/process/MucfProcess.cc new file mode 100644 index 0000000000..cb2947e0a6 --- /dev/null +++ b/src/celeritas/mucf/process/MucfProcess.cc @@ -0,0 +1,67 @@ +//------------------------------- -*- C++ -*- -------------------------------// +// Copyright Celeritas contributors: see top-level COPYRIGHT file for details +// SPDX-License-Identifier: (Apache-2.0 OR MIT) +//---------------------------------------------------------------------------// +//! \file celeritas/mucf/process/MucfProcess.cc +//---------------------------------------------------------------------------// +#include "MucfProcess.hh" + +#include + +#include "celeritas/mat/MaterialParams.hh" +#include "celeritas/mucf/model/DTMixMucfModel.hh" +#include "celeritas/phys/Model.hh" +#include "celeritas/phys/ParticleParams.hh" + +namespace celeritas +{ +//---------------------------------------------------------------------------// +/*! + * Construct from host data. + */ +MucfProcess::MucfProcess(SPConstParticles particles, SPConstMaterials materials) + : particles_(particles), materials_(materials) +{ + //! \todo Fix ImportProcessClass + CELER_EXPECT(particles_); + CELER_EXPECT(materials_); +} + +//---------------------------------------------------------------------------// +/*! + * Construct the models associated with this process. + */ +auto MucfProcess::build_models(ActionIdIter start_id) const -> VecModel +{ + return {std::make_shared( + *start_id++, *particles_, *materials_)}; +} + +//---------------------------------------------------------------------------// +/*! + * Get the interaction cross sections for the given energy range. + */ +auto MucfProcess::macro_xs(Applicability) const -> XsGrid +{ + return {}; +} + +//---------------------------------------------------------------------------// +/*! + * Get the energy loss for the given energy range. + */ +auto MucfProcess::energy_loss(Applicability) const -> EnergyLossGrid +{ + return {}; +} +//---------------------------------------------------------------------------// +/*! + * Name of the process. + */ +std::string_view MucfProcess::label() const +{ + return "Muon-catalyzed fusion"; +} + +//---------------------------------------------------------------------------// +} // namespace celeritas diff --git a/src/celeritas/mucf/process/MucfProcess.hh b/src/celeritas/mucf/process/MucfProcess.hh new file mode 100644 index 0000000000..23949c2acc --- /dev/null +++ b/src/celeritas/mucf/process/MucfProcess.hh @@ -0,0 +1,61 @@ +//------------------------------- -*- C++ -*- -------------------------------// +// Copyright Celeritas contributors: see top-level COPYRIGHT file for details +// SPDX-License-Identifier: (Apache-2.0 OR MIT) +//---------------------------------------------------------------------------// +//! \file celeritas/mucf/process/MucfProcess.hh +//---------------------------------------------------------------------------// +#pragma once + +#include "celeritas/Types.hh" +#include "celeritas/phys/Process.hh" + +namespace celeritas +{ +class MaterialParams; +class ParticleParams; + +//---------------------------------------------------------------------------// +/*! + * Muon-catalyzed fusion of muonic dd, dt, or tt molecules. + * + * \note This is an at-rest process. + */ +class MucfProcess final : public Process +{ + public: + //!@{ + //! \name Type aliases + using SPConstParticles = std::shared_ptr; + using SPConstMaterials = std::shared_ptr; + //!@} + + public: + // Construct from shared and imported data + explicit MucfProcess(SPConstParticles particles, + SPConstMaterials materials); + + // Construct the models associated with this process + VecModel build_models(ActionIdIter start_id) const final; + + // Get the interaction cross sections for the given energy range + XsGrid macro_xs(Applicability) const final; + + // Get the energy loss for the given energy range + EnergyLossGrid energy_loss(Applicability) const final; + + //! Whether the integral method can be used to sample interaction length + bool supports_integral_xs() const final { return false; } + + //! Whether the process applies when the particle is stopped + bool applies_at_rest() const final { return true; } + + // Name of the process + std::string_view label() const final; + + private: + SPConstParticles particles_; + SPConstMaterials materials_; +}; + +//---------------------------------------------------------------------------// +} // namespace celeritas diff --git a/src/celeritas/phys/PDGNumber.hh b/src/celeritas/phys/PDGNumber.hh index 44755468c1..15baf0edbf 100644 --- a/src/celeritas/phys/PDGNumber.hh +++ b/src/celeritas/phys/PDGNumber.hh @@ -146,6 +146,10 @@ CELER_DEFINE_PDGNUMBER(anti_he3, -1000020030) CELER_DEFINE_PDGNUMBER(alpha, 1000020040) CELER_DEFINE_PDGNUMBER(anti_alpha, -1000020040) +// Muonic nuclei +CELER_DEFINE_PDGNUMBER(muonic_deuteron, 2000010020) +CELER_DEFINE_PDGNUMBER(muonic_triton, 2000010030) + //!@} #undef CELER_DEFINE_PDGNUMBER diff --git a/test/celeritas/data/four-steel-slabs.root-dump.json b/test/celeritas/data/four-steel-slabs.root-dump.json index 22948401fa..cba069f76c 100644 --- a/test/celeritas/data/four-steel-slabs.root-dump.json +++ b/test/celeritas/data/four-steel-slabs.root-dump.json @@ -569,6 +569,25 @@ }, "mucf_physics" : { "_typename" : "celeritas::inp::MucfPhysics", + "scalars" : { + "_typename" : "celeritas::inp::MucfScalars", + "protium" : { + "_typename" : "celeritas::Quantity", + "value_" : 0 + }, + "deuterium" : { + "_typename" : "celeritas::Quantity", + "value_" : 0 + }, + "tritium" : { + "_typename" : "celeritas::Quantity", + "value_" : 0 + }, + "liquid_hydrogen_density" : { + "_typename" : "celeritas::Quantity", + "value_" : 0 + } + }, "muon_energy_cdf" : { "_typename" : "celeritas::inp::Grid", "x" : [], diff --git a/test/celeritas/ext/GeantImporter.test.cc b/test/celeritas/ext/GeantImporter.test.cc index 68022e0a32..a42578192a 100644 --- a/test/celeritas/ext/GeantImporter.test.cc +++ b/test/celeritas/ext/GeantImporter.test.cc @@ -2132,15 +2132,27 @@ TEST_F(OpticalSurfaces, surfaces) } //---------------------------------------------------------------------------// -TEST_F(MucfBox, run) +TEST_F(MucfBox, static_data) { auto const& mucf = this->imported_data().mucf_physics; EXPECT_TRUE(mucf); - static double const expected_muon_energy_cdf_y[] = {1, 1}; - EXPECT_EQ(2, mucf.muon_energy_cdf.x.size()); - EXPECT_VEC_EQ(expected_muon_energy_cdf_y, mucf.muon_energy_cdf.y); + auto average = [](inp::Grid::VecDbl const& data) -> double { + double sum{0}; + for (auto y : data) + { + sum += y; + } + return sum / data.size(); + }; + + static size_t const expected_muon_energy_cdf_size = 21; + EXPECT_EQ(expected_muon_energy_cdf_size, mucf.muon_energy_cdf.x.size()); + EXPECT_EQ(expected_muon_energy_cdf_size, mucf.muon_energy_cdf.y.size()); + EXPECT_SOFT_EQ(0.55157437567861023, average(mucf.muon_energy_cdf.x)); + EXPECT_SOFT_EQ(11.250286274435437, average(mucf.muon_energy_cdf.y)); + // Dummy data auto const& cycle_f0 = mucf.cycle_rates[0]; static double const expected_cycle_rate_f0_y[] = {2, 2}; EXPECT_TRUE(cycle_f0);