[WIP] Python improvements and output_dir implementation

examples/make_input_file_and_run.py

@@ -255,16 +255,22 @@
     'geometry_input': geometry_0D
 }
 
+# Test rustbca_py and 'output_dir' option
 input_data['options']['name'] = 'rustbca_input_file'
+if not os.path.exists(r'./outputs'):
+    os.makedirs('outputs')
+input_data['options']['output_dir'] = 'outputs'
 rustbca_py(input_data, mode)
-s = np.genfromtxt('rustbca_input_filesputtered.output', delimiter=',')
+s = np.genfromtxt('outputs/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])
 
+# reset these options before running from command line
 input_data['options']['name'] = 'input_file'
+input_data['options']['output_dir'] = '.'
 
 # Attempt to cleanup line endings
 input_string = dumps(input_data).replace('\r', '')

src/consts.rs

@@ -44,4 +44,6 @@ pub const TABLE_SIZE: usize = Z_MAX*(Z_MAX + 1)/2;
 /// Gauss-Legendre Points
 pub const GL_X: [f64; 5] = [0./2. + 1./2., -0.538469/2. + 1./2., 0.538469/2. + 1./2., -0.90618/2. + 1./2., 0.90618/2. + 1./2.];
 /// Gauss-Legendre Weights
-pub const GL_W: [f64; 5] = [0.568889/2., 0.478629/2., 0.478629/2., 0.236927/2., 0.236927/2.];
+pub const GL_W: [f64; 5] = [0.568889/2., 0.478629/2., 0.478629/2., 0.236927/2., 0.236927/2.];
+/// Enforced maximum number density
+pub const MAX_DENSITY: f64 = 1e32;

src/geometry.rs

@@ -65,8 +65,9 @@ impl Geometry for Mesh0D {
 
         let densities: Vec<f64> = input.densities.iter().map(|element| element/(length_unit).powi(3)).collect();
         assert!(!densities.is_empty(), "Input Error: density list empty.");
-
+        
         let total_density: f64 = densities.iter().sum();
+        assert!(total_density < MAX_DENSITY, "Input Error: total density {}/m^3 exceeds realistic values; check values or units.", total_density);
 
         let energy_barrier_thickness = 1./total_density.cbrt()/SQRTPI*2.;
 
@@ -154,6 +155,7 @@ impl Geometry for Mesh1D {
         let densities: Vec<Vec<f64>> = geometry_input.densities
             .iter()
             .map( |row| row.iter().map(|element| element/(length_unit).powi(3)).collect() ).collect();
+
 
         //Assert all layer density lists are equal length
         assert!(
@@ -170,6 +172,7 @@ impl Geometry for Mesh1D {
             layer_bottom += layer_thickness*length_unit;
 
             let total_density: f64 = densities.iter().sum();
+            assert!(total_density < MAX_DENSITY, "Input Error: total density {}/m^3 exceeds realistic values; check values or units.", total_density);
             let concentrations: Vec<f64> = densities.iter().map(|&density| density/total_density).collect::<Vec<f64>>();
 
             layers.push(Layer1D::new(layer_top, layer_bottom, densities, concentrations, ck));
@@ -308,6 +311,7 @@ impl Geometry for HomogeneousMesh2D {
         let densities: Vec<f64> = input.densities.iter().map(|element| element/(length_unit).powi(3)).collect();
 
         let total_density: f64 = densities.iter().sum();
+        assert!(total_density < MAX_DENSITY, "Input Error: total density {}/m^3 exceeds realistic values; check values or units.", total_density);
 
         let energy_barrier_thickness = 1./total_density.cbrt()/SQRTPI*2.;
 
@@ -476,6 +480,7 @@ impl Geometry for Mesh2D {
                 y3*length_unit,
             );
             let total_density: f64 = densities.iter().sum();
+            assert!(total_density < MAX_DENSITY, "Input Error: total density {}/m^3 exceeds realistic values; check values or units.", total_density);
             let concentrations: Vec<f64> = densities.iter().map(|&density| density/total_density).collect::<Vec<f64>>();
 
             cells.push(Cell2D::new(coordinate_set_converted, densities, concentrations, ck));

src/input.rs

@@ -212,6 +212,10 @@ fn default_rootfinder() -> Vec<Vec<Rootfinder>> {
     vec![vec![Rootfinder::DEFAULTNEWTON]]
 }
 
+fn default_output_dir() -> String {
+    "./".to_string()
+}
+
 /// Rustbca's internal representation of the simulation-level options.
 #[cfg(not(feature = "distributions"))]
 #[derive(Deserialize, Clone)]
@@ -250,7 +254,9 @@ pub struct Options {
     #[serde(default = "default_false")]
     pub track_energy_losses: bool,
     #[serde(default = "default_seed")]
-    pub seed: i32
+    pub seed: i32,
+    #[serde(default = "default_output_dir")]
+    pub output_dir: String,
 }
 
 #[cfg(not(feature = "distributions"))]
@@ -275,6 +281,7 @@ impl Options {
             track_displacements: false,
             track_energy_losses: false,
             seed: default_seed(),
+            output_dir: default_output_dir(),
         }
     }
 }
@@ -331,7 +338,9 @@ pub struct Options {
     pub y_num: usize,
     pub z_num: usize,
     #[serde(default = "default_seed")]
-    pub seed: i32
+    pub seed: i32,
+    #[serde(default = "default_output_dir")]
+    pub output_dir: String,
 }
 
 #[cfg(feature = "distributions")]
@@ -370,7 +379,8 @@ impl Options {
             x_num: 0,
             y_num: 0,
             z_num: 0,
-            seed: default_seed()
+            seed: default_seed(),
+            output_dir: default_output_dir(),
         }
     }
 }

src/lib.rs

@@ -1196,19 +1196,19 @@ pub fn compound_bca_list_tracked_py<'py>(energies: Vec<f64>, ux: Vec<f64>, uy: V
 ///Returns:
 ///    vx, vy, vz (float): final x, y, and z velocity in m/s. When ion implants in material, vx, vy, and vz will all be zero.
 #[pyfunction]
-pub fn reflect_single_ion_py<'py>(ion: &Bound<'py, PyDict>, target: &Bound<'py, PyDict>, vx: f64, vy: f64, vz: f64) -> (f64, f64, f64){
+pub fn reflect_single_ion_py<'py>(ion: &Bound<'py, PyDict>, target: &Bound<'py, PyDict>, vx: f64, vy: f64, vz: f64) -> PyResult<(f64, f64, f64)> {
 
-    let Z1: f64 = ion.get_item("Z").unwrap().expect("Error: Cannot get key 'Z' from ion dict.").extract().unwrap();
-    let m1: f64 = ion.get_item("m").unwrap().expect("Error: Cannot get key 'm' from ion dict.").extract().unwrap();
-    let Es1: f64 = ion.get_item("Es").unwrap().expect("Error: Cannot get key 'Es' from ion dict.").extract().unwrap();
-    let Ec1: f64 = ion.get_item("Ec").unwrap().expect("Error: Cannot get key 'Ec' from ion dict.").extract().unwrap();
-
-    let Z2: f64 = target.get_item("Z").unwrap().expect("Error: Cannot get key 'Z' from target dict.").extract().unwrap();
-    let m2: f64 = target.get_item("m").unwrap().expect("Error: Cannot get key 'm' from target dict.").extract().unwrap();
-    let Es2: f64 = target.get_item("Es").unwrap().expect("Error: Cannot get key 'Es' from target dict.").extract().unwrap();
-    let Ec2: f64 = target.get_item("Ec").unwrap().expect("Error: Cannot get key 'Ec' from target dict.").extract().unwrap();
-    let Eb2: f64 = target.get_item("Eb").unwrap().expect("Error: Cannot get key 'Eb' from target dict.").extract().unwrap();
-    let n2: f64 = target.get_item("n").unwrap().expect("Error: Cannot get key 'n' from target dict.").extract().unwrap();
+    let Z1: f64 = get_value_from_dict!(ion, "Z")?;
+    let m1: f64 = get_value_from_dict!(ion, "m")?;
+    let Es1: f64 = get_value_from_dict!(ion, "Es")?;
+    let Ec1: f64 = get_value_from_dict!(ion, "Ec")?;
+
+    let Z2: f64 = get_value_from_dict!(target, "Z")?;
+    let m2: f64 = get_value_from_dict!(target, "m")?;
+    let Es2: f64 = get_value_from_dict!(target, "Es")?;
+    let Ec2: f64 = get_value_from_dict!(target, "Ec")?;
+    let Eb2: f64 = get_value_from_dict!(target, "Eb")?;
+    let n2: f64 = get_value_from_dict!(target, "n")?;
 
     assert!(vx > 0.0, "Input error: vx must be greater than zero for incident particles to hit surface at x=0.");
 
@@ -1271,9 +1271,9 @@ pub fn reflect_single_ion_py<'py>(ion: &Bound<'py, PyDict>, target: &Bound<'py,
     let vz2 = output[0].dir.z*reflected_velocity;
 
     if output[0].E > 0.0 && output[0].dir.x < 0.0 && output[0].left && output[0].incident {
-        (vx2, vy2, vz2)
+        Ok((vx2, vy2, vz2))
     } else {
-        (0.0, 0.0, 0.0)
+        Ok((0.0, 0.0, 0.0))
     }
 }
 

src/material.rs

@@ -229,17 +229,6 @@ impl <T: Geometry> Material<T> {
         }
     }
 
-    ///The minimum cutoff energy of all species that make up the material.
-    pub fn minimum_cutoff_energy(&self) -> f64 {
-        let mut min_Ec = self.Ec.iter().sum::<f64>();
-        for Ec in self.Ec.iter() {
-            if min_Ec > *Ec {
-                min_Ec = *Ec;
-            }
-        }
-        min_Ec
-    }
-
     ///Choose the parameters of a target atom as a concentration-weighted random draw from the species in the triangle that contains or is nearest to (x, y).
     pub fn choose(&self, x: f64, y: f64, z: f64, rng: &mut ChaCha8Rng) -> (usize, f64, f64, f64, f64, f64, usize) {
         let random_number: f64 = rng.random::<f64>();

src/output.rs

@@ -1,5 +1,25 @@
 use super::*;
 use std::fs::File;
+use std::path::Path;
+#[cfg(feature = "distributions")]
+extern crate ndarray;
+#[cfg(feature = "distributions")]
+use ndarray::prelude::*;
+
+macro_rules! open_output_file {
+    ($output_path:expr, $name:expr, $options:expr) => {{
+    //Open output files for streaming output
+    let file_path = $output_path.join(format!("{}{}{}", $options.name, $name, ".output"));
+    let file = OpenOptions::new()
+        .write(true)
+        .create(true)
+        .truncate(true)
+        .open(file_path)
+        .context(format!("Could not open {} output file.", $options.name,))
+        .unwrap();
+    BufWriter::with_capacity($options.write_buffer_size, file)
+    }}
+}
 
 #[derive(Clone, Debug)]
 pub struct OutputUnits {
@@ -21,12 +41,6 @@ pub fn energy_angle_from_particle(particle: &particle::Particle, units: &OutputU
     (energy, angle)
 }
 
-#[cfg(feature = "distributions")]
-extern crate ndarray;
-
-#[cfg(feature = "distributions")]
-use ndarray::prelude::*;
-
 /// Distribution tracker for tracking EADs and implantation distributions
 #[derive(Serialize)]
 #[cfg(feature = "distributions")]
@@ -67,17 +81,22 @@ impl Distributions {
     /// Write distributions to toml
     pub fn print(&self, options: &Options) {
 
+        let output_path = Path::new(&options.output_dir);
+        assert!(output_path.try_exists().unwrap());
+        assert!(output_path.is_dir());
+
+        let file_path = output_path.join(format!("{}{}{}", options.name, "distributions", ".toml"));
         let distribution_output_file = OpenOptions::new()
             .write(true)
             .create(true)
             .truncate(true)
-            .open(format!("{}{}", options.name, "distributions.toml"))
+            .open(file_path)
             .context("Could not open distributions output file.")
             .unwrap();
-        let mut distribution_file_stream = BufWriter::with_capacity(8000, distribution_output_file);
+        let mut distribution_file_stream = BufWriter::with_capacity(options.write_buffer_size, distribution_output_file);
         let toml = toml::to_string(&self).unwrap();
+
         writeln!(distribution_file_stream, "{}", toml).unwrap();
-
     }
 
     /// Updates distributions with a single particle
@@ -203,14 +222,15 @@ pub struct SummaryPerSpecies {
 impl SummaryPerSpecies {
     pub fn new(options: &Options) -> SummaryPerSpecies {
 
-        let summary_output_file = OpenOptions::new()
-            .write(true)
-            .create(true)
-            .truncate(true)
-            .open(format!("{}{}", options.name, "summary.output"))
-            .context("Could not open output file.")
-            .unwrap();
-        let writer = BufWriter::with_capacity(8000, summary_output_file);
+        let output_path = Path::new(&options.output_dir);
+        assert!(output_path.try_exists().unwrap());
+        assert!(output_path.is_dir());
+
+        let writer = open_output_file!(
+            output_path,
+            "summary",
+            options
+        );
 
         SummaryPerSpecies {
             m: vec![],
@@ -259,69 +279,53 @@ impl SummaryPerSpecies {
 
 /// Open list output files for streaming write
 pub fn open_output_lists(options: &Options) -> OutputListStreams {
-    //Open output files for streaming output
-    let reflected_file = OpenOptions::new()
-        .write(true)
-        .create(true)
-        .truncate(true)
-        .open(format!("{}{}", options.name, "reflected.output"))
-        .context("Could not open output file.")
-        .unwrap();
-    let reflected_file_stream = BufWriter::with_capacity(options.write_buffer_size, reflected_file);
 
-    let sputtered_file = OpenOptions::new()
-        .write(true)
-        .create(true)
-        .truncate(true)
-        .open(format!("{}{}", options.name, "sputtered.output"))
-        .context("Could not open output file.")
-        .unwrap();
-    let sputtered_file_stream = BufWriter::with_capacity(options.write_buffer_size, sputtered_file);
+    let output_path = Path::new(&options.output_dir);
+    assert!(output_path.try_exists().unwrap());
+    assert!(output_path.is_dir());
 
-    let deposited_file = OpenOptions::new()
-        .write(true)
-        .create(true)
-        .truncate(true)
-        .open(format!("{}{}", options.name, "deposited.output"))
-        .context("Could not open output file.")
-        .unwrap();
-    let deposited_file_stream = BufWriter::with_capacity(options.write_buffer_size, deposited_file);
-
-    let trajectory_file = OpenOptions::new()
-        .write(true)
-        .create(true)
-        .truncate(true)
-        .open(format!("{}{}", options.name, "trajectories.output"))
-        .context("Could not open output file.")
-        .unwrap();
-    let trajectory_file_stream = BufWriter::with_capacity(options.write_buffer_size, trajectory_file);
-
-    let trajectory_data = OpenOptions::new()
-        .write(true)
-        .create(true)
-        .truncate(true)
-        .open(format!("{}{}", options.name, "trajectory_data.output"))
-        .context("Could not open output file.")
-        .unwrap();
-    let trajectory_data_stream = BufWriter::with_capacity(options.write_buffer_size, trajectory_data);
-
-    let displacements_file = OpenOptions::new()
-        .write(true)
-        .create(true)
-        .truncate(true)
-        .open(format!("{}{}", options.name, "displacements.output"))
-        .context("Could not open output file.")
-        .unwrap();
-    let displacements_file_stream = BufWriter::with_capacity(options.write_buffer_size, displacements_file);
-
-    let energy_loss_file = OpenOptions::new()
-        .write(true)
-        .create(true)
-        .truncate(true)
-        .open(format!("{}{}", options.name, "energy_loss.output"))
-        .context("Could not open output file.")
-        .unwrap();
-    let energy_loss_file_stream = BufWriter::with_capacity(options.write_buffer_size, energy_loss_file);
+    //Open output files for streaming output
+    let reflected_file_stream = open_output_file!(
+        output_path,
+        "reflected",
+        options
+    );
+
+    let sputtered_file_stream = open_output_file!(
+        output_path,
+        "sputtered",
+        options
+    );
+
+    let deposited_file_stream = open_output_file!(
+        output_path,
+        "deposited",
+        options
+    );
+
+    let trajectory_file_stream = open_output_file!(
+        output_path,
+        "trajectories",
+        options
+    );
+
+    let trajectory_data_stream = open_output_file!(
+        output_path,
+        "trajectory_data",
+        options
+    );
+
+    let displacements_file_stream = open_output_file!(
+        output_path,
+        "displacements",
+        options
+    );
+
+    let energy_loss_file_stream = open_output_file!(
+        output_path,
+        "energy_loss",
+        options
+    );
 
     OutputListStreams {
         reflected_file_stream,