From b810a83fd696c0421defea495baf420091e0bd1b Mon Sep 17 00:00:00 2001 From: Jon Drobny <37962344+drobnyjt@users.noreply.github.com> Date: Mon, 10 Aug 2026 17:51:42 -0700 Subject: [PATCH 1/4] Update README.md --- README.md | 13 ++++++------- 1 file changed, 6 insertions(+), 7 deletions(-) diff --git a/README.md b/README.md index d196d01c..c25a3930 100644 --- a/README.md +++ b/README.md @@ -134,6 +134,7 @@ The following features are implemented in `RustBCA`: * nonlocal (Lindhard-Scharff), * and equipartition * Biersack-Varelas interpolation is also included for electronic stopping up to ~1 GeV/nucleon. Note that high energy physics beyond electronic stopping are not included, and that Biersack-Varelas may not be as accurate as other methods. +* A modified Biersack-Varelas form which includes an interpolation parameter, `c_i`, that can be used to better match available stopping data for the Bragg peak. * Biersack-Haggmark treatment of high-energy free-flight paths between collisions can be included to greatly speed up high-energy simulations (i.e., by neglecting very small angle scattering). * A wide range of interaction potentials are provided, including: * the Kr-C, ZBL, Lenz-Jensen, and Moliere universal, screened-Coulomb potentials. @@ -157,7 +158,8 @@ The following features are implemented in `RustBCA`: * full trajectory tracking for both the incident ions and target atoms, * and many other parameters such as position of origin of sputtered particles and energy loss along trajectories. * Optionally, the code can produce energy-angle and implantation distributions when built with the `--features distributions` flag and disable space-intensive particle list output with `--features no_list_output`. -* Library functions for modeling ion reflection, implantation, and sputtering in C++/C, Python, and Fortran codes. +* Library functions for running RustBCA in-memory in C++/C, Python, and Fortran codes. +* A seeded PRNG such that RustBCA results are exactly reproducible on any machine for any number of threads ## Installation @@ -170,8 +172,6 @@ cargo build --release will add an executable at `target/release/`. -[HDF5] for particle list input has been tested on Windows, but version 1.10.6 must be used. - #### Manual Dependences * [rustup], the [Rust] toolchain (includes `cargo`, the [Rust] package manager, `rustc`, the [Rust] compiler, and more). @@ -182,7 +182,6 @@ will add an executable at `target/release/`. #### Optional Dependencies -* [HDF5] libraries * For manipulating input files and running associated scripts, the following are suggested: * [Python] 3.6+ * [Python] libraries: `numpy`, `matplotlib`, `toml`, `shapely`, and `scipy`. @@ -216,9 +215,9 @@ git clone https://github.com/lcpp-org/RustBCA cd RustBCA cargo build --release ``` -8. (Optional) Build `RustBCA` with optional dependencies, `hdf5` and/or `rcpr`: +8. (Optional) Build `RustBCA` with optional dependencies such as `rcpr`: ```bash -cargo build --release --features cpr_rootfinder,hdf5 +cargo build --release --features cpr_rootfinder ``` 9. `input.toml` is the input file - see the [Input File](https://github.com/lcpp-org/RustBCA/wiki/Standalone-Code:-Input-File) page for more information 10. Run the required tests using: @@ -284,7 +283,7 @@ Additionally, `RustBCA` accepts an input file type (one of: `0D`, `1D`, `2D`, `T ```bash ./RustBCA 0D /path/to/input.toml ``` -**Warning: RustBCA defaults to the 2D triangular mesh input mode.** For more details, see [Input Files](https://github.com/lcpp-org/RustBCA/wiki/Standalone-Code:-Input-File). +**Note: RustBCA defaults to the 2D triangular mesh input mode.** For more details, see [Input Files](https://github.com/lcpp-org/RustBCA/wiki/Standalone-Code:-Input-File). Also have a look at the examples on the [Wiki] to see some examples of RustBCA input files. [BCA]: https://en.wikipedia.org/wiki/Binary_collision_approximation From ac6faa98d6254f1fbed2f83fdbfedeb8208f3afc Mon Sep 17 00:00:00 2001 From: Jon Drobny Date: Mon, 10 Aug 2026 19:01:17 -0700 Subject: [PATCH 2/4] Updated curated citations in readme. --- README.md | 7 ++++++- 1 file changed, 6 insertions(+), 1 deletion(-) diff --git a/README.md b/README.md index c25a3930..f0168c01 100644 --- a/README.md +++ b/README.md @@ -23,7 +23,11 @@ Journal of Open Source Software by clicking the badge below: [![DOI](https://joss.theoj.org/papers/10.21105/joss.03298/status.svg)](https://doi.org/10.21105/joss.03298) -Selected citations of RustBCA as of 5/24/23: +Selected citations of RustBCA as of 8/10/26: +* [Validation of hybrid-PIC Simulations for Advanced Beam-Driven FRC Modeling](https://iopscience.iop.org/article/10.1088/1741-4326/ae96c0/meta), R. E. Groenewald, et al. (2026) +* [(Preprint): Comparative qualification of advanced plasma-facing materials for fusion pilot plants through public- and private-sector experiments in DIII-D](https://arxiv.org/abs/2607.23400), F. Effenberg, et al. (2026) +* [Understanding carbon sourcing and transport originating from the helicon antenna surfaces during high-power helicon discharge in DIII-D Tokamak](https://iopscience.iop.org/article/10.1088/1741-4326/ae9322/meta), A. Kumarm et al. (2026) +* [Molecular dynamics simulations of reflection and sputtering behavior of boron under deuterium ion irradiation](https://www.sciencedirect.com/science/article/pii/S0022311526001327?casa_token=axebkNZe94gAAAAA:WG47RTEOepmq2gZnvPnpFWGfsLHOY1gibRY-Ed0wLXl6zEyUkQPFZnOGrRJhtR-3SVkh57p62Z4) H. Schamis, et al. (2026) * [Simulation of liquid lithium divertor geometry using SOLPS-ITER](https://doi.org/10.1109/TPS.2022.3166402), JD Lore et al. (2022) * [Characterizing W sources in the all-W wall, all-RF WEST tokamak environment](https://doi.org/10.1088/1361-6587/ac8acc), CC Klepper et al. (2022) * [hPIC2: A hardware-accelerated, hybrid particle-in-cell code for dynamic plasma-material interactions](https://doi.org/10.1016/j.cpc.2022.108569), LT Meredith et al. (2023) @@ -31,6 +35,7 @@ Selected citations of RustBCA as of 5/24/23: * [Modeling the effect of nitrogen recycling on the erosion and leakage of tungsten impurities from the SAS-VW divertor in DIII-D during nitrogen gas injection](https://doi.org/10.1016/j.nme.2022.101254), MS Parsons et al. (2023) * [Enabling attractive-repulsive potentials in binary-collision-approximation monte-carlo codes for ion-surface interactions](https://doi.org/10.1088/2053-1591/ad1262), J Drobny and D Curreli (2023) * [Multi-physics modeling of tungsten collector probe samples during the WEST C4 He campaign](https://doi.org.10.1088/1741-4326/ad6c5b), A. Lasa et al. (2024) +* [Integrated modeling of RF-induced tungsten erosion at ICRH antenna structures in the WEST tokamak*](https://iopscience.iop.org/article/10.1088/1741-4326/ade455/meta), A. Kumar et al., (2025) ## Getting started From 00fe5d00a48826f78eb3181b58536a9d2f3d7507 Mon Sep 17 00:00:00 2001 From: Jon Drobny Date: Mon, 10 Aug 2026 22:19:06 -0700 Subject: [PATCH 3/4] Converted sputtering_yield to a rayon sum; reflection coefficient will be harder to incorporate the residue-preserving sum... --- src/lib.rs | 47 +++++++++++++++++++++-------------------------- 1 file changed, 21 insertions(+), 26 deletions(-) diff --git a/src/lib.rs b/src/lib.rs index 39703eab..4419c7a9 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -1811,38 +1811,33 @@ pub fn sputtering_yield<'py>(ion: &Bound<'py, PyDict>, target: &Bound<'py, PyDic let x = -m.geometry.energy_barrier_thickness; - let num_sputtered = Mutex::new(0); - let seed: u64 = get_seed().map_err(|error| PyValueError::new_err(""))?; - (0..num_samples as u64).into_par_iter() - .for_each_init( + Ok((0..num_samples as u64).into_par_iter() + .map_init( || ChaCha8Rng::seed_from_u64(seed), |rng, index| { - let p = particle::Particle::default_incident( - m1, - Z1, - energy, - Ec1, - Es1, - x, - ux, - uy, - uz - ); - - rng.set_stream(index); - let output = bca::single_ion_bca(p, &m, &options, rng); + let p = particle::Particle::default_incident( + m1, + Z1, + energy, + Ec1, + Es1, + x, + ux, + uy, + uz + ); + - for particle in output { - if particle.E > 0.0 && particle.dir.x < 0.0 && particle.left && (!particle.incident) { - let mut num_sputtered = num_sputtered.lock().unwrap(); - *num_sputtered += 1; - } + rng.set_stream(index); + let output = bca::single_ion_bca(p, &m, &options, rng); + + output.iter() + .filter(|particle| particle.E > 0.0 && particle.dir.x < 0.0 && particle.left && !particle.incident) + .count() } - }); - let num_sputtered = *num_sputtered.lock().unwrap(); - Ok(num_sputtered as f64 / num_samples as f64) + ).sum::() as f64 / num_samples as f64) } #[cfg(feature = "python")] From d36b2f2c39e57350107fb366d8737120c5972101 Mon Sep 17 00:00:00 2001 From: Jon Drobny Date: Tue, 11 Aug 2026 09:19:42 -0700 Subject: [PATCH 4/4] Removed mutexes in lib.rs; rayon sums preserve order so no longer need TwoSum --- src/lib.rs | 137 ++++++++++++++++++++--------------------------------- 1 file changed, 52 insertions(+), 85 deletions(-) diff --git a/src/lib.rs b/src/lib.rs index 4419c7a9..02b7b2ae 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -30,9 +30,6 @@ use std::os::raw::c_int; //standard slice use std::slice; -//Mutex for multithreading in ergonomic Python library functions -#[cfg(feature = "python")] -use std::sync::Mutex; //itertools use itertools::{izip}; @@ -1905,54 +1902,38 @@ pub fn reflection_coefficient<'py>(ion: &Bound<'py, PyDict>, target: &Bound<'py, let x = -m.geometry.energy_barrier_thickness; - let num_reflected = Mutex::new(0); - let energy_reflected = Mutex::new(0.0); - let residue = Mutex::new(0.0); - let seed: u64 = get_seed().map_err(|error| PyValueError::new_err(""))?; - (0..num_samples as u64).into_par_iter() - .for_each_init( + let (num_reflected, energy_reflected): (Vec, Vec) = (0..num_samples as u64).into_par_iter() + .map_init( || ChaCha8Rng::seed_from_u64(seed), |rng, index| { - let p = particle::Particle::default_incident( - m1, - Z1, - energy, - Ec1, - Es1, - x, - ux, - uy, - uz - ); - - rng.set_stream(index); - let output = bca::single_ion_bca(p, &m, &options, rng); - - for particle in output { - if particle.E > 0.0 && particle.dir.x < 0.0 && particle.left && particle.incident { - let mut num_reflected = num_reflected.lock().unwrap(); - *num_reflected += 1; - - let mut energy_reflected = energy_reflected.lock().unwrap(); - - let residue_part; - - // Use Moller-Knuth TwoSum to preserve deterministic fp reduce - (*energy_reflected, residue_part) = moller_knuth_two_sum(*energy_reflected, particle.E); + let p = particle::Particle::default_incident( + m1, + Z1, + energy, + Ec1, + Es1, + x, + ux, + uy, + uz + ); + + rng.set_stream(index); + let output = bca::single_ion_bca(p, &m, &options, rng); - let mut residue = residue.lock().unwrap(); - *residue = *residue + residue_part; + let mut count = 0; + let mut energy = 0.0; + for particle in output.into_iter().filter(|particle| particle.E > 0.0 && particle.dir.x < 0.0 && particle.left && particle.incident) { + count += 1; + energy += particle.E; } - }; - }); - if let (Ok(num_reflected), Ok(energy_reflected), Ok(residue)) = (num_reflected.lock(), energy_reflected.lock(), residue.lock()) { - return Ok((*num_reflected as f64 / num_samples as f64, (*energy_reflected + *residue) / EV / energy / num_samples as f64)) - } else { - return Err(PyValueError::new_err("Check input values.")) - } - + (count, energy) + } + ).collect(); + // map_init preserves order - that means we don't have to use TwoSum for reproducibility + Ok((num_reflected.into_iter().sum::() as f64 / num_samples as f64, energy_reflected.iter().sum::() / EV / energy / num_samples as f64)) } fn get_seed() -> Result { @@ -2061,52 +2042,38 @@ pub fn compound_reflection_coefficient<'py>(ion: &Bound<'py, PyDict>, targets: V let x = -m.geometry.energy_barrier_thickness; - let num_reflected = Mutex::new(0); - let energy_reflected = Mutex::new(0.0); - let residue = Mutex::new(0.0); - let seed: u64 = get_seed().map_err(|error| PyValueError::new_err(""))?; - (0..num_samples as u64).into_par_iter() - .for_each_init( + let (num_reflected, energy_reflected): (Vec, Vec) = (0..num_samples as u64).into_par_iter() + .map_init( || ChaCha8Rng::seed_from_u64(seed), |rng, index| { + let p = particle::Particle::default_incident( + m1, + Z1, + energy, + Ec1, + Es1, + x, + ux, + uy, + uz + ); + + rng.set_stream(index); + let output = bca::single_ion_bca(p, &m, &options, rng); - let p = particle::Particle::default_incident( - m1, - Z1, - energy, - Ec1, - Es1, - x, - ux, - uy, - uz - ); - - rng.set_stream(index); - let output = bca::single_ion_bca(p, &m, &options, rng); - - for particle in output { - if particle.E > 0.0 && particle.dir.x < 0.0 && particle.left && particle.incident { - let mut num_reflected = num_reflected.lock().unwrap(); - *num_reflected += 1; - let mut energy_reflected = energy_reflected.lock().unwrap(); - - let residue_part; - - // Use Moller-Knuth TwoSum to preserve deterministic fp reduce - (*energy_reflected, residue_part) = moller_knuth_two_sum(*energy_reflected, particle.E); - - let mut residue = residue.lock().unwrap(); - *residue = *residue + residue_part; + let mut count = 0; + let mut energy = 0.0; + for particle in output.into_iter().filter(|particle| particle.E > 0.0 && particle.dir.x < 0.0 && particle.left && particle.incident) { + count += 1; + energy += particle.E; } - } - }); - let num_reflected = *num_reflected.lock().unwrap(); - let energy_reflected = *energy_reflected.lock().unwrap(); - let residue = *residue.lock().unwrap(); - Ok((num_reflected as f64 / num_samples as f64, (energy_reflected + residue) / EV / energy / num_samples as f64)) + (count, energy) + } + ).collect(); + // map_init preserves order - that means we don't have to use TwoSum for reproducibility + Ok((num_reflected.into_iter().sum::() as f64 / num_samples as f64, energy_reflected.iter().sum::() / EV / energy / num_samples as f64)) } /// Moller-Knuth TwoSum Floating-Point Adder with Residual (FPAR)