diff --git a/Cargo.toml b/Cargo.toml index 6b2c374..6933a13 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -29,6 +29,7 @@ rcpr = { git = "https://github.com/drobnyjt/rcpr", optional = true} ndarray = {version = "0.17.2", features = ["serde"], optional = true} parry3d-f64 = {optional = true, version="0.2.0"} pyo3 = {version = "0.29.0", optional=true} +pythonize = {version = "0.29.0", optional=true} [dev-dependencies] float-cmp = "0.10.0" @@ -45,5 +46,5 @@ cpr_rootfinder = ["rcpr"] distributions = ["ndarray"] no_list_output = [] parry3d = ["parry3d-f64"] -python = ["pyo3"] +python = ["pyo3", "pythonize"] extended_max_z = [] \ No newline at end of file diff --git a/examples/make_input_file_and_run.py b/examples/make_input_file_and_run.py index b48d376..37e3a16 100644 --- a/examples/make_input_file_and_run.py +++ b/examples/make_input_file_and_run.py @@ -255,6 +255,17 @@ 'geometry_input': geometry_0D } +input_data['options']['name'] = 'rustbca_input_file' +rustbca_py(input_data, mode) +s = np.genfromtxt('rustbca_input_filesputtered.output', delimiter=',') + +arrays = rustbca_local_py(input_data, mode) +sputtered = arrays['sputtered'] + +np.testing.assert_approx_equal(s[0, 2], np.array(arrays['energy'])[sputtered][0]) + +input_data['options']['name'] = 'input_file' + # Attempt to cleanup line endings input_string = dumps(input_data).replace('\r', '') with open('examples/input_file.toml', 'w') as input_file: @@ -265,6 +276,9 @@ # Read output files - ensure arrays are at least 2D for indexing sputtered = np.atleast_2d(np.genfromtxt('input_filesputtered.output', delimiter=',')) + +np.testing.assert_approx_equal(sputtered[0, 2], s[0, 2]) + reflected = np.atleast_2d(np.genfromtxt('input_filereflected.output', delimiter=',')) implanted = np.atleast_2d(np.genfromtxt('input_filedeposited.output', delimiter=',')) diff --git a/examples/test_different_options.py b/examples/test_different_options.py index 587d6ce..f3ca629 100644 --- a/examples/test_different_options.py +++ b/examples/test_different_options.py @@ -51,6 +51,8 @@ a TOML file using tomlkit. It runs the input file with cargo run --release and reads the output files. + +It also runs with rustbca_py and ensures the output files are identical. ''' def run_test( @@ -275,6 +277,19 @@ def run_test( reflected = np.atleast_2d(np.genfromtxt(f'input_file_{index}reflected.output', delimiter=',')) implanted = np.atleast_2d(np.genfromtxt(f'input_file_{index}deposited.output', delimiter=',')) + input_data['options']['name'] = f'python_{index}' + if run_sim: + rustbca_py(input_data, geometry_mode=mode) + + # Read output files - ensure arrays are at least 2D for indexing + sputtered_py = np.atleast_2d(np.genfromtxt(f'python_{index}sputtered.output', delimiter=',')) + reflected_py = np.atleast_2d(np.genfromtxt(f'python_{index}reflected.output', delimiter=',')) + implanted_py = np.atleast_2d(np.genfromtxt(f'python_{index}deposited.output', delimiter=',')) + + np.testing.assert_allclose(sputtered, sputtered_py) + np.testing.assert_allclose(reflected, reflected_py) + np.testing.assert_allclose(implanted, implanted_py) + return sputtered, reflected, implanted num_bins = 75 diff --git a/examples/test_scattering_integrals.py b/examples/test_scattering_integrals.py index 7e5d9c4..2a43da1 100644 --- a/examples/test_scattering_integrals.py +++ b/examples/test_scattering_integrals.py @@ -10,7 +10,7 @@ from formulas import * energies = np.logspace(0, 4, 4) -impact_parameters = np.logspace(-3, 3, 100) +impact_parameters = np.logspace(-3, 2, 100) ion = helium target = boron @@ -24,23 +24,28 @@ linestyles = ['-', '--', ':', '-.'] -for linestyle, energy in zip(linestyles, energies): - gm = np.zeros(100) - gl = np.zeros(100) - mw = np.zeros(100) - magic = np.zeros(100) - for index, p in enumerate(impact_parameters): - gm[index], gl[index], mw[index], magic[index] = scattering_integrals(Za, Zb, Ma, Mb, energy, p) - plt.semilogx(impact_parameters, gm, label=f'Gauss-Mehler, E={np.round(energy/1000)} keV', linestyle=linestyle) - plt.semilogx(impact_parameters, gl, label=f'Gauss-Legendre, E={np.round(energy/1000)} keV', linestyle=linestyle) - plt.semilogx(impact_parameters, mw, label=f'Mendenhall-Weller, E={np.round(energy/1000)} keV', linestyle=linestyle) - plt.semilogx(impact_parameters, magic, label=f'MAGIC, E={np.round(energy/1000)} keV', linestyle=linestyle) - plt.gca().set_prop_cycle(None) - - np.testing.assert_allclose(gm, gl, atol=5e-3) # 0.5% seems reasonable? max is ~0.3% - np.testing.assert_allclose(gm, mw, atol=5e-3) - np.testing.assert_allclose(mw, gl, atol=5e-3) - plt.legend() - plt.xlabel('p [A]') - plt.ylabel('theta [rad]') +for potential in ["KR_C", "MOLIERE", "ZBL"]: + plt.figure() + plt.title(f'Scattering Angles for {potential}') + for linestyle, energy in zip(linestyles, energies): + gm = np.zeros(100) + gl = np.zeros(100) + mw = np.zeros(100) + magic = np.zeros(100) + for index, p in enumerate(impact_parameters): + gm[index], gl[index], mw[index], magic[index] = scattering_integrals(Za, Zb, Ma, Mb, energy, p, interaction_potential=potential) + + plt.semilogx(impact_parameters, gm, label=f'Gauss-Mehler, E={np.round(energy/1000, 3)} keV', linestyle=linestyle) + plt.semilogx(impact_parameters, gl, label=f'Gauss-Legendre, E={np.round(energy/1000, 3)} keV', linestyle=linestyle) + plt.semilogx(impact_parameters, mw, label=f'Mendenhall-Weller, E={np.round(energy/1000, 3)} keV', linestyle=linestyle) + plt.semilogx(impact_parameters, magic, label=f'MAGIC, E={np.round(energy/1000, 3)} keV', linestyle=linestyle) + plt.gca().set_prop_cycle(None) + + np.testing.assert_allclose(gm, gl, atol=5e-3) # 0.5% seems reasonable? max is ~0.3% + np.testing.assert_allclose(gm, mw, atol=5e-3) + np.testing.assert_allclose(mw, gl, atol=5e-3) + plt.legend() + plt.xlabel('p [A]') + plt.ylabel('theta [rad]') + if show_plots: plt.show() \ No newline at end of file diff --git a/src/input.rs b/src/input.rs index 02665b1..e6c607f 100644 --- a/src/input.rs +++ b/src/input.rs @@ -381,9 +381,7 @@ impl Options { } } -pub fn input(input_file: String) -> (Vec, material::Material, Options, OutputUnits) -where ::InputFileFormat: Deserialize<'static> + 'static { - +pub fn read_input_file(input_file: String) -> ::InputFileFormat { //Read input file, convert to string, and open with toml let mut input_toml = String::new(); let mut file = OpenOptions::new() @@ -394,7 +392,19 @@ where ::InputFileFormat: Deserialize<'static> + 'static { .unwrap_or_else(|_| panic!("Input errror: could not open input file {}.", &input_file)); file.read_to_string(&mut input_toml).context("Could not convert TOML file to string.").unwrap(); - let input: ::InputFileFormat = InputFile::new(&input_toml); + InputFile::new(&input_toml) +} + +pub fn input(input_file: String) -> (Vec, material::Material, Options, OutputUnits) +where ::InputFileFormat: Deserialize<'static> + 'static { + + let input: ::InputFileFormat = read_input_file::(input_file); + + process_input_file(input) + +} + +pub fn process_input_file(input: ::InputFileFormat) -> (Vec, material::Material, Options, OutputUnits) { //Unpack toml information into structs let options = (*input.get_options()).clone(); diff --git a/src/lib.rs b/src/lib.rs index 9d24cae..67ebc4a 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -9,10 +9,10 @@ use std::alloc::{dealloc, Layout}; use std::mem::align_of; //Parallelization - currently only used in python library functions +//#[cfg(feature = "python")] +//use rayon::ThreadPoolBuilder; #[cfg(feature = "python")] -use rayon::prelude::*; -#[cfg(feature = "python")] -use rayon::*; +use rayon::iter::{IndexedParallelIterator, ParallelExtend, IntoParallelIterator, ParallelIterator}; //Error handling crate use anyhow::{Result, Context, anyhow}; @@ -53,6 +53,10 @@ use std::f64::consts::SQRT_2; use pyo3::prelude::*; #[cfg(feature = "python")] use pyo3::types::*; +#[cfg(feature = "python")] +use pythonize::*; +#[cfg(feature = "python")] +use pyo3::exceptions::{PyValueError, PyRuntimeError}; //Load internal modules pub mod material; @@ -68,9 +72,7 @@ pub mod consts; pub mod structs; pub mod sphere; pub mod math; - -#[cfg(feature = "parry3d")] -pub mod parry; +pub mod physics; pub use crate::enums::*; pub use crate::consts::*; @@ -81,6 +83,10 @@ pub use crate::geometry::{Geometry, GeometryElement, Mesh0D, Mesh1D, Mesh2D}; pub use crate::sphere::{Sphere, SphereInput, InputSphere}; pub use crate::math::*; pub use crate::material::*; +pub use crate::physics::*; + +#[cfg(feature = "parry3d")] +pub mod parry; #[cfg(feature = "parry3d")] pub use crate::parry::{ParryBall, ParryBallInput, InputParryBall, ParryTriMesh, ParryTriMeshInput, InputParryTriMesh}; @@ -139,6 +145,12 @@ mod libRustBCA { #[pymodule_export] use super::scattering_integrals; + + #[pymodule_export] + use super::rustbca_py; + + #[pymodule_export] + use super::rustbca_local_py; } #[derive(Debug)] @@ -2131,7 +2143,6 @@ pub fn compound_reflection_coefficient<'py>(ion: &Bound<'py, PyDict>, targets: V let mut residue = residue.lock().unwrap(); *residue = *residue + residue_part; - } } }); @@ -2161,20 +2172,118 @@ fn moller_knuth_two_sum(a: f64, b: f64) -> (f64, f64) { let r = delta_a + delta_b; (s, r) } + #[cfg(feature = "python")] #[pyfunction] -#[pyo3(signature = (Za, Zb, Ma, Mb, E0, p, n_gl_points=100))] -fn scattering_integrals(Za: f64, Zb: f64, Ma: f64, Mb: f64, E0: f64, p: f64, n_gl_points: usize) -> (f64, f64, f64, f64) { +#[pyo3(signature = (Za, Zb, Ma, Mb, E0, p, n_gl_points=100, interaction_potential="KR_C"))] +fn scattering_integrals(Za: f64, Zb: f64, Ma: f64, Mb: f64, E0: f64, p: f64, n_gl_points: usize, interaction_potential: &str) -> PyResult<(f64, f64, f64, f64)> { let E0 = E0*EV; let p = p*ANGSTROM; - let x0_newton = bca::newton_rootfinder(Za, Zb, Ma, Mb, E0, p, InteractionPotential::KR_C, 1000, 1E-12).unwrap(); + let potential = match interaction_potential { + "KR_C" => InteractionPotential::KR_C, + "LENZ_JENSEN" => InteractionPotential::LENZ_JENSEN, + "MOLIERE" => InteractionPotential::MOLIERE, + "ZBL" => InteractionPotential::ZBL, + _ => return Err(PyValueError::new_err(format!("Unimplemented interaction potential {}; try 'KR_C'", interaction_potential))) + }; + + let x0_newton = bca::newton_rootfinder(Za, Zb, Ma, Mb, E0, p, potential, 1000, 1E-12).map_err( + |error| PyRuntimeError::new_err(format!("Rootfinder failed to find distance of closest approach; check input values.")) + )?; //Compute center of mass deflection angle with each algorithm - let theta_gm = bca::gauss_mehler(Za, Zb, Ma, Mb, E0, p, x0_newton, InteractionPotential::KR_C, n_gl_points); - let theta_gl = bca::gauss_legendre(Za, Zb, Ma, Mb, E0, p, x0_newton, InteractionPotential::KR_C); - let theta_mw = bca::mendenhall_weller(Za, Zb, Ma, Mb, E0, p, x0_newton, InteractionPotential::KR_C); - let theta_magic = bca::magic(Za, Zb, Ma, Mb, E0, p, x0_newton, InteractionPotential::KR_C); + let theta_gm = bca::gauss_mehler(Za, Zb, Ma, Mb, E0, p, x0_newton, potential, n_gl_points); + let theta_gl = bca::gauss_legendre(Za, Zb, Ma, Mb, E0, p, x0_newton, potential); + let theta_mw = bca::mendenhall_weller(Za, Zb, Ma, Mb, E0, p, x0_newton, potential); + let theta_magic = bca::magic(Za, Zb, Ma, Mb, E0, p, x0_newton, potential); + + Ok((theta_gm, theta_gl, theta_mw, theta_magic)) +} +#[cfg(feature = "python")] +macro_rules! geometry_typed_loops { + ($geometry_type:ty, $input:expr, $python:expr) => { + { + let input: <$geometry_type as geometry::Geometry>::InputFileFormat = depythonize(&$input).unwrap(); + let (particle_input_array, material, options, output_units) = input::process_input_file(input); + let pool = rayon::ThreadPoolBuilder::new().num_threads(options.num_threads).build().unwrap(); + pool.install( || + physics::physics_loop::<$geometry_type>(particle_input_array, material, options, output_units) + ); + Ok(()) + } + } +} + +#[cfg(feature = "python")] +#[pyfunction] +#[pyo3(signature=(input, geometry_mode="1D"))] +fn rustbca_py<'py>(python: Python<'py>, input: &Bound<'py, PyDict>, geometry_mode: &str) -> PyResult<()> { + match geometry_mode { + "0D" => geometry_typed_loops!(Mesh0D, input, python), + "1D" => geometry_typed_loops!(Mesh1D, input, python), + "2D" => geometry_typed_loops!(Mesh2D, input, python), + "HOMOGENEOUS2D" => geometry_typed_loops!(Mesh2D, input, python), + "SPHERE" => geometry_typed_loops!(Sphere, input, python), + #[cfg(feature="parry3d")] + "BALL" => geometry_typed_loops!(ParryBall, input, python), + #[cfg(feature="parry3d")] + "TRIMESH" => geometry_typed_loops!(ParryTriMesh, input, python), + _ => Err(PyValueError::new_err(format!("Input Error: Unimplemented geometry mode {}; try '1D'", geometry_mode))) + } +} - (theta_gm, theta_gl, theta_mw, theta_magic) +/* +Notes on macros - this is the first I have written, so I'm taking notes here as I go. +macro_rules! makes a macro - here, the macro is called geometry_types_silent_loops +macros pattern match an argument and replace it with anything you want +I want it to take a tuple of a string (e.g., "1D") and a type (e.g., Mesh1D) +and plop those into corresponding match arms. +The first line tells the macro to expect an argument with that pattern. +arguments are $:. Designators: +block +expr is used for expressions +ident is used for variable/function names +item +literal is used for literal constants +pat (pattern) +path +stmt (statement) +tt (token tree) +ty (type) +vis (visibility qualifier) +*/ +#[cfg(feature = "python")] +macro_rules! geometry_typed_silent_loops { + ($geometry_type:ty, $input:expr, $python:expr) => { + { + let input: <$geometry_type as geometry::Geometry>::InputFileFormat = depythonize(&$input).unwrap(); + let (particle_input_array, material, options, output_units) = input::process_input_file(input); + let pool = rayon::ThreadPoolBuilder::new().num_threads(options.num_threads).build().unwrap(); + let finished_particles = pool.install( || + physics::silent_physics_loop::<$geometry_type>(particle_input_array, material, options, output_units.clone()) + ); + let finished_particles_container = physics::process_finished_particles_to_arrays(finished_particles, output_units); + Ok(pythonize($python, &finished_particles_container)?) + } + } +} + +#[cfg(feature = "python")] +#[pyfunction] +#[pyo3(signature=(input, geometry_mode="1D"))] +fn rustbca_local_py<'py>(python: Python<'py>, input: &Bound<'py, PyDict>, geometry_mode: &str) -> PyResult> { + + match geometry_mode { + "0D" => geometry_typed_silent_loops!(Mesh0D, input, python), + "1D" => geometry_typed_silent_loops!(Mesh1D, input, python), + "2D" => geometry_typed_silent_loops!(Mesh2D, input, python), + "HOMOGENEOUS2D" => geometry_typed_silent_loops!(Mesh2D, input, python), + "SPHERE" => geometry_typed_silent_loops!(Sphere, input, python), + #[cfg(feature="parry3d")] + "BALL" => geometry_typed_silent_loops!(ParryBall, input, python), + #[cfg(feature="parry3d")] + "TRIMESH" => geometry_typed_silent_loops!(ParryTriMesh, input, python), + _ => Err(PyValueError::new_err(format!("Input Error: Unimplemented geometry mode {}; try '1D'", geometry_mode))) + } } \ No newline at end of file diff --git a/src/main.rs b/src/main.rs index 8b02d4b..0920812 100644 --- a/src/main.rs +++ b/src/main.rs @@ -5,9 +5,6 @@ use std::{env, fmt}; use std::mem::discriminant; -//Progress bar crate - works with rayon -use indicatif::{ProgressBar, ProgressStyle}; - //Error handling crate use anyhow::{Result, Context, anyhow}; @@ -19,8 +16,8 @@ use serde::*; use hdf5::*; //Parallelization -use rayon::prelude::*; -use rayon::*; +//use rayon::prelude::*; +//use rayon::ThreadPoolBuilder; //I/O use std::fs::OpenOptions; @@ -69,6 +66,20 @@ pub use crate::math::duff_orthonormal_basis; #[cfg(feature = "parry3d")] pub use crate::parry::{ParryBall, ParryBallInput, InputParryBall, ParryTriMesh, ParryTriMeshInput, InputParryTriMesh}; + +macro_rules! main_loop { + ($geometry_type:ident, $input_file:expr) => { + { + let (particle_input_array, material, options, output_units) = input::input::<$geometry_type>($input_file); + //Initialize threads with rayon + println!("Processing {} ions...", particle_input_array.len()); + println!("Initializing with {} threads...", options.num_threads); + let _ = rayon::ThreadPoolBuilder::new().num_threads(options.num_threads).build_global(); + physics_loop::<$geometry_type>(particle_input_array, material, options, output_units); + } + } +} + fn main() { let args: Vec = env::args().collect(); @@ -91,36 +102,16 @@ fn main() { _ => panic!("Too many command line arguments. RustBCA accepts 0 (use 'input.toml') 1 () or 2 ( )"), }; - match geometry_type { - GeometryType::MESH0D => { - let (particle_input_array, material, options, output_units) = input::input::(input_file); - physics_loop::(particle_input_array, material, options, output_units); - }, - GeometryType::MESH1D => { - let (particle_input_array, material, options, output_units) = input::input::(input_file); - physics_loop::(particle_input_array, material, options, output_units); - }, - GeometryType::MESH2D => { - let (particle_input_array, material, options, output_units) = input::input::(input_file); - physics_loop::(particle_input_array, material, options, output_units); - }, - GeometryType::SPHERE => { - let (particle_input_array, material, options, output_units) = input::input::(input_file); - physics_loop::(particle_input_array, material, options, output_units); - }, + // This invokes the above macro that expands into the physics loop invocation for each type + match geometry_type { + GeometryType::MESH0D => main_loop!(Mesh0D, input_file), + GeometryType::MESH1D => main_loop!(Mesh1D, input_file), + GeometryType::MESH2D => main_loop!(Mesh2D, input_file), + GeometryType::SPHERE => main_loop!(Sphere, input_file), #[cfg(feature = "parry3d")] - GeometryType::BALL => { - let (particle_input_array, material, options, output_units) = input::input::(input_file); - physics_loop::(particle_input_array, material, options, output_units); - } + GeometryType::BALL => main_loop!(ParryBall, input_file), #[cfg(feature = "parry3d")] - GeometryType::TRIMESH => { - let (particle_input_array, material, options, output_units) = input::input::(input_file); - physics_loop::(particle_input_array, material, options, output_units); - } - GeometryType::HOMOGENEOUS2D => { - let (particle_input_array, material, options, output_units) = input::input::(input_file); - physics_loop::(particle_input_array, material, options, output_units); - } + GeometryType::TRIMESH => main_loop!(ParryTriMesh, input_file), + GeometryType::HOMOGENEOUS2D => main_loop!(HomogeneousMesh2D, input_file), } } diff --git a/src/output.rs b/src/output.rs index ad3d7a2..463c0c5 100644 --- a/src/output.rs +++ b/src/output.rs @@ -1,6 +1,7 @@ use super::*; use std::fs::File; +#[derive(Clone, Debug)] pub struct OutputUnits { pub length_unit: f64, pub energy_unit: f64, diff --git a/src/physics.rs b/src/physics.rs index 933d358..177b911 100644 --- a/src/physics.rs +++ b/src/physics.rs @@ -1,8 +1,87 @@ use super::*; +use indicatif::{ProgressBar, ProgressStyle}; +use rayon::iter::{IndexedParallelIterator, ParallelExtend, IntoParallelIterator, ParallelIterator}; + +pub fn silent_physics_loop(particle_input_array: Vec, material: material::Material, options: Options, output_units: OutputUnits) -> Vec { + + let mut finished_particles: Vec = Vec::new(); + + finished_particles.par_extend( + particle_input_array.into_par_iter() + .enumerate() + .map_init( + || if options.seed < 0 { + ChaCha8Rng::seed_from_u64(rand::random()) + } else { + ChaCha8Rng::seed_from_u64(u64::try_from(options.seed).expect("Value Error: seed not u64.")) + }, + | rng, (particle_index, particle_input)| { + rng.set_stream((particle_index) as u64); + bca::single_ion_bca(particle::Particle::from_input(particle_input, &options), &material, &options, rng) + }).flatten() + ); + finished_particles +} -pub fn physics_loop(particle_input_array: Vec, material: material::Material, options: Options, output_units: OutputUnits) { +#[derive(Clone, Debug, Serialize, Deserialize)] +pub struct FinishedParticlesContainer { + sputtered: Vec, + implanted: Vec, + atomic_number: Vec, + mass: Vec, + energy: Vec, + x: Vec, + y: Vec, + z: Vec, + ux: Vec, + uy: Vec, + uz: Vec, +} +impl FinishedParticlesContainer { + pub fn new() -> FinishedParticlesContainer { + FinishedParticlesContainer { + sputtered: Vec::new(), + implanted: Vec::new(), + atomic_number: Vec::new(), + mass: Vec::new(), + energy: Vec::new(), + x: Vec::new(), + y: Vec::new(), + z: Vec::new(), + ux: Vec::new(), + uy: Vec::new(), + uz: Vec::new(), + } + } +} - println!("Processing {} ions...", particle_input_array.len()); +pub fn process_finished_particles_to_arrays(finished_particles: Vec, output_units: OutputUnits) -> FinishedParticlesContainer { + let mut finished_particles_container = FinishedParticlesContainer::new(); + for particle in finished_particles { + + let implanted = particle.incident & !particle.left; + let sputtered = !particle.incident & particle.left; + let reflected = particle.incident & particle.left; + + if implanted | sputtered | reflected { + finished_particles_container.sputtered.push(sputtered); + finished_particles_container.implanted.push(implanted); + + finished_particles_container.atomic_number.push(particle.Z as usize); + finished_particles_container.mass.push(particle.m/output_units.mass_unit); + finished_particles_container.energy.push(particle.E/output_units.energy_unit); + finished_particles_container.x.push(particle.pos.x/output_units.length_unit); + finished_particles_container.y.push(particle.pos.y/output_units.length_unit); + finished_particles_container.z.push(particle.pos.z/output_units.length_unit); + finished_particles_container.ux.push(particle.dir.x); + finished_particles_container.uy.push(particle.dir.y); + finished_particles_container.uz.push(particle.dir.z); + } + } + finished_particles_container +} + +pub fn physics_loop(particle_input_array: Vec, material: material::Material, options: Options, output_units: OutputUnits) { let total_count: u64 = particle_input_array.len() as u64; assert!(total_count/options.num_chunks > 0, "Input error: chunk size == 0 - reduce num_chunks or increase particle count."); @@ -15,10 +94,6 @@ pub fn physics_loop(particle_input_array: Vec