Skip to content

Remove mutexes from lib.rs #333

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
drobnyjt merged 5 commits into from
Aug 11, 2026
Merged

Remove mutexes from lib.rs #333

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
5 commits
Select commit Hold shift + click to select a range
673bd3f
Merge pull request #328 from lcpp-org/dev
drobnyjt Aug 11, 2026
b810a83
Update README.md
drobnyjt Aug 11, 2026
ac6faa9
Updated curated citations in readme.
drobnyjt Aug 11, 2026
00fe5d0
Converted sputtering_yield to a rayon sum; reflection coefficient wil…
drobnyjt Aug 11, 2026
d36b2f2
Removed mutexes in lib.rs; rayon sums preserve order so no longer nee…
drobnyjt Aug 11, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
20 changes: 12 additions & 8 deletions README.md
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -23,14 +23,19 @@ 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)
* [Global sensitivity analysis of a coupled multiphysics model to predict surface evolution in fusion plasma–surface interactions](https://doi.org/10.1016/j.commatsci.2023.112229), P. Robbe et al. (2023)
* [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

Expand Down Expand Up @@ -134,6 +139,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.
Expand All @@ -157,7 +163,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

Expand All @@ -170,8 +177,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).
Expand All @@ -182,7 +187,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`.
Expand Down Expand Up @@ -216,9 +220,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:
Expand Down Expand Up @@ -284,7 +288,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
Expand Down
188 changes: 75 additions & 113 deletions src/lib.rs
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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};
Expand Down Expand Up @@ -1811,38 +1808,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::<usize>() as f64 / num_samples as f64)
}

#[cfg(feature = "python")]
Expand Down Expand Up @@ -1910,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<usize>, Vec<f64>) = (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 mut residue = residue.lock().unwrap();
*residue = *residue + residue_part;
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 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::<usize>() as f64 / num_samples as f64, energy_reflected.iter().sum::<f64>() / EV / energy / num_samples as f64))
}

fn get_seed() -> Result<u64> {
Expand Down Expand Up @@ -2066,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<usize>, Vec<f64>) = (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 mut residue = residue.lock().unwrap();
*residue = *residue + residue_part;
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 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::<usize>() as f64 / num_samples as f64, energy_reflected.iter().sum::<f64>() / EV / energy / num_samples as f64))
}

/// Moller-Knuth TwoSum Floating-Point Adder with Residual (FPAR)
Expand Down
Loading