Light ion transport #3752
Closed
Light ion transport #3752
+4,236
−39
Conversation
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
- Extend ParticleType enum with proton, deuteron, triton, helion, alpha - Add IncidentChargedParticle class for loading HDF5 cross section data - Implement log-log interpolation for cross section lookup - Add particle_type_from_string and particle_type_to_string conversions - Update energy_min/energy_max arrays from size 4 to 9 - Add C++ unit test for charged particle data loading - Add Python script to generate HDF5 from ENDF files Co-Authored-By: Claude Opus 4.5 <noreply@anthropic.com>
Add proton_transport, deuteron_transport, triton_transport, helion_transport, and alpha_transport settings to both C++ and Python. - Add settings to settings.h and settings.cpp - Add XML parsing in read_settings_xml - Reset settings in finalize.cpp - Recognize charged particle sources in source.cpp - Add Python properties and XML serialization in settings.py Co-Authored-By: Claude Opus 4.5 <noreply@anthropic.com>
- Add ejectile particle type to IncidentChargedParticle class - Store ejectile types in HDF5 as integer indices - Add charged_particle_ index vector to Material class - Implement sample_charged_particle_reaction with proper two-body kinematics - Handle particle type changes when ejectile differs from projectile - Add mass constants for charged particles (proton, deuteron, triton, helion, alpha) - Update particle_type_to_str to include charged particle types - Add ejectile type verification in unit tests Co-Authored-By: Claude Opus 4.5 <noreply@anthropic.com>
…ions Replace classical kinematics with relativistic formulas from ENDF-6 Appendix E for charged particle reaction sampling. Key changes: - Add TwoBodyKinematicsResult struct for ejectile and residual data - Implement relativistic_two_body_kinematics() following ENDF Eq. E.4-E.27 - Add identify_particle_from_mass() to determine residual particle type - Update sample_charged_particle_reaction() to create both ejectile and residual particles when transport is enabled - Add D-T fusion unit test verifying ~14.1 MeV neutron and ~3.5 MeV alpha The relativistic formulas correctly handle mass-energy conversion via Q-value, giving residual mass m_R = m_t + (m_i - m_e) - Q. Energy conservation is verified to ~10^-12 keV precision. Co-Authored-By: Claude Opus 4.5 <noreply@anthropic.com>
This commit adds the continuous slowing down approximation (CSDA) for charged particle transport: - Add stopping power models to Material (constant and SRIM-based) - Add charged_particle_path setting for cross section data location - Add charged_particle_delta setting for fractional energy loss per step - Implement CSDA transport loop in event_advance_charged_particle() - Auto-load charged particle data files based on transport settings - Add energy cutoff validation for charged particles - Fix nuclide sampling to recalculate XS at collision energy - Add verbose output for charged particle transport (verbosity >= 10) - Add regression test for D-T fusion The CSDA loop samples collision distance at each substep, properly handling the interplay between energy loss and collision probability. Co-Authored-By: Claude Opus 4.5 <noreply@anthropic.com>
When neutrons undergo elastic scattering with light nuclei (H1, H2, H3, He3, He4) that have charged particle transport enabled, the recoil nucleus is now created as a secondary particle. This is important for fusion applications where 14 MeV neutrons can impart significant kinetic energy to deuterons and tritons, which may then undergo fusion reactions themselves. Changes: - Add nuclide_to_particle_type() helper to map nuclide names to particle types - Add is_charged_particle_transport_enabled() helper function - Modify elastic_scatter() to calculate recoil energy using momentum conservation and create secondary particles above the energy cutoff - Add bounds checking for charged_particle_ vector access in CSDA loop - Add verbose output for elastic recoil creation (verbosity >= 10) Co-Authored-By: Claude Opus 4.5 <noreply@anthropic.com>
- Add --targets option to process only specific target nuclides - Add --subdirs flag to optionally create particle-type subdirectories - Default to flat output directory (all HDF5 files in one folder) - Add docstring documentation for process_particle_files function Co-Authored-By: Claude Opus 4.5 <noreply@anthropic.com>
Add analytic Bethe-Bloch stopping power calculation using PDG tabulated mean excitation energies for elements Z=1-92. The implementation includes: - PDG mean excitation energy table from adndt.pdf - Bragg additivity rule for effective I in compound materials - Relativistic kinematics (beta, gamma, T_max) - Low-energy clamping when formula breaks down Validated against NIST PSTAR data for protons in beryllium: - 10 MeV: 0.3% difference - 1 MeV: 4% difference - Lower energies show expected deviations due to missing shell corrections Also adds electron_density() method to Material class, default Nuclide constructor for testing, and extends verbosity range to level 11. Co-Authored-By: Claude Opus 4.5 <noreply@anthropic.com>
This commit adds support for energy-dependent product yields in charged particle reactions, which is essential for correctly modeling reactions like MT=5 (lumped reactions) that can produce multiple particle types with yields that vary with incident energy. Key changes: - Add NJOY ACER templates for charged particle ACE file generation - Add from_ace() method to IncidentChargedParticle for loading ACE data - Extend HDF5 format to include reaction products with yields - Add products_ member to IncidentChargedParticle C++ class - Modify sample_charged_particle_reaction to use product yields - Add ParticleType::unknown for residual nuclei that aren't transported - Extend str_to_particle_type to handle GNDS names (H1, H2, He3, etc.) The proton-Be9 test case now correctly produces ~0.0034 neutron current (previously ~0.0012 due to missing neutrons from MT=5). Co-Authored-By: Claude Opus 4.5 <noreply@anthropic.com>
Two bugs were causing segfaults when sampling secondary particle distributions for charged particle reactions: 1. UncorrelatedAngleEnergy::sample() unconditionally dereferenced energy_ even when null. This occurs when HDF5 data contains only angular distributions without energy data (e.g., deuteron_H3.h5 for D+T fusion). Added null check with fallback to incident energy. 2. NBodyPhaseSpace::sample() only handled n=3,4,5 bodies and threw a misleading error for other values. The deuteron data uses n=2 for some products. Added: - Case for n=2 (uniform energy distribution) - General case for arbitrary n using ENDF-6 Eq. 6.21 - Helper function sample_gamma_half() for Gamma(0.5,1) sampling Co-Authored-By: Claude Opus 4.5 <noreply@anthropic.com>
Add a new global tally that tracks total neutron production from secondary particle creation. This is useful for fusion and other applications where neutron yield is a key metric. Changes: - Add NEUTRON_PRODUCTION to GlobalTally enum - Increment N_GLOBAL_TALLIES from 4 to 5 - Add neutron_production_tally_ member to particle data - Score neutron weight in create_secondary() when type is neutron - Accumulate to global tally in event_death() with atomic pragma - Report "Neutron Production" in results output with uncertainties The tally is normalized per source particle, similar to leakage fraction. Co-Authored-By: Claude Opus 4.5 <noreply@anthropic.com>
In eigenvalue mode, all secondary neutrons are now placed in the fission bank instead of the secondary bank, allowing neutron multiplication chains (e.g., from (n,2n) reactions) to propagate across generations. This enables k-eigenvalue calculations for systems without fissile material. The implementation applies the same 1/keff population control as fission: - Score k-tracklength based on actual neutron weight - Sample number of sites from expected value (weight/keff) - Create sites with unit weight Co-Authored-By: Claude Opus 4.5 <noreply@anthropic.com>
Add the Li-Petrasso stopping power formula for charged particles in a plasma, which accounts for both electron and ion stopping with proper treatment of the Chandrasekhar function and Coulomb logarithm. The implementation includes: - Chandrasekhar psi function and its derivative - Electron stopping with velocity-dependent psi(x) - Ion stopping with the G function for equal-mass scattering - Debye length and quantum minimum impact parameter - Temperature-dependent plasma physics This enables simulation of charged particle transport in hot plasmas where traditional Bethe-Bloch cold matter stopping is not applicable. Reference: C. K. Li and R. D. Petrasso, Phys. Rev. Lett. 70, 3059 (1993) Co-Authored-By: Claude Opus 4.5 <noreply@anthropic.com>
Only bank secondary neutrons when the source particle is not a neutron (e.g., from charged particle fusion reactions). Neutrons produced by other neutrons (via (n,2n) etc.) go to the secondary bank as usual. This restores the original eigenvalue mode semantics where the fission bank tracks multiplying reactions, not the full neutron chain. Co-Authored-By: Claude Opus 4.5 <noreply@anthropic.com>
Instead of assuming hydrogen plasma (n_e = rho/m_p), use the material's actual electron density from electron_density() method. This correctly handles arbitrary material compositions. The electron_density() returns e/b-cm, which is converted to cm^-3 by multiplying by 1e24 for use in the CGS-based Li-Petrasso formula. Co-Authored-By: Claude Opus 4.5 <noreply@anthropic.com>
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
1 participant
Add this suggestion to a batch that can be applied as a single commit.
This suggestion is invalid because no changes were made to the code.
Suggestions cannot be applied while the pull request is closed.
Suggestions cannot be applied while viewing a subset of changes.
Only one suggestion per line can be applied in a batch.
Add this suggestion to a batch that can be applied as a single commit.
Applying suggestions on deleted lines is not supported.
You must change the existing code in this line in order to create a valid suggestion.
Outdated suggestions cannot be applied.
This suggestion has been applied or marked resolved.
Suggestions cannot be applied from pending reviews.
Suggestions cannot be applied on multi-line comments.
Suggestions cannot be applied while the pull request is queued to merge.
Suggestion cannot be applied right now. Please check back later.
Description
This is a relatively complete implementation of light ion transport I wanted to use to explore the possibility of a chain reaction of DT fusion -> 14.1 MeV neutrons -> elastic scattering producing low MeV D's and T's -> more fusion. Turns out the multiplication factor for this both in infinite medium solids and plasmas is like... < 0.01. So it could not be used to build anything workable. It was inspired by this paper.
I am simply opening and then closing this PR since I don't really have the time to make this clean and merge-able, but would like to put the code out there for others to potentially cherry pick from and use. It is by no means complete or full-featured, but in a basic validation test case for a range of proton energies incident on beryllium seems to give answers in solid agreement with MCNP. As a disclaimer this is vibe-coded but with considerable user intervention and a bit of validation.
An example of a shortcoming is that Q<0 reactions might not be handled correctly at the moment (ctrl-f for TODO in the diff here). Another potential shortcoming is that I'm not sure if the angular distributions are correct at all. There are also some extraneous changes in here like an interpretation of k-eigenvalue with DT fusion looking like fission that would be applicable in that aforementioned hypothesized chain reaction, and a global tally for total neutron production per source particle.
The MCNP data was digitized from this recent paper.
Checklist