From 7f8f01850f009aebf39111e1630170c1cd188e4f Mon Sep 17 00:00:00 2001 From: Jon Drobny Date: Sat, 27 Feb 2021 09:45:55 -0600 Subject: [PATCH 01/14] First commit of geometry rewrite. Mesh and MeshElement traits have been created. First steps towards geometry agnosticism. --- src/bca.rs | 34 ++++++++-------- src/main.rs | 1 + src/material.rs | 105 ++++++++++++++++++++++++++---------------------- src/mesh.rs | 51 +++++++++++++++-------- src/tests.rs | 28 ++++++------- 5 files changed, 124 insertions(+), 95 deletions(-) diff --git a/src/bca.rs b/src/bca.rs index 0fea7b7..624e6dc 100644 --- a/src/bca.rs +++ b/src/bca.rs @@ -96,7 +96,7 @@ pub fn single_ion_bca(particle: particle::Particle, material: &material::Materia &binary_collision_geometry, &options); //If recoil location is inside, proceed with binary collision loop - if material.inside(particle_2.pos.x, particle_2.pos.y) & material.inside_energy_barrier(particle_1.pos.x, particle_1.pos.y) { + if material.inside(particle_2.pos.x, particle_2.pos.y, particle_2.pos.z) & material.inside_energy_barrier(particle_1.pos.x, particle_1.pos.y, particle_1.pos.z) { //Determine scattering angle from binary collision let binary_collision_result = bca::calculate_binary_collision(&particle_1, @@ -113,7 +113,7 @@ pub fn single_ion_bca(particle: particle::Particle, material: &material::Materia } //Energy transfer to recoil - particle_2.E = binary_collision_result.recoil_energy - material.average_bulk_binding_energy(particle_2.pos.x, particle_2.pos.y); + particle_2.E = binary_collision_result.recoil_energy - material.average_bulk_binding_energy(particle_2.pos.x, particle_2.pos.y, particle_2.pos.z); particle_2.energy_origin = particle_2.E; //Accumulate asymptotic deflections for primary particle @@ -128,7 +128,7 @@ pub fn single_ion_bca(particle: particle::Particle, material: &material::Materia particle::rotate_particle(&mut particle_2, -binary_collision_result.psi_recoil, binary_collision_geometry.phi_azimuthal); - + particle_2.dir_old.x = particle_2.dir.x; particle_2.dir_old.y = particle_2.dir.y; particle_2.dir_old.z = particle_2.dir.z; @@ -148,13 +148,14 @@ pub fn single_ion_bca(particle: particle::Particle, material: &material::Materia let Ma: f64 = particle_1.m; let Mb: f64 = particle_2.m; - let n = material.total_number_density(particle_2.pos.x, particle_2.pos.y); + let n = material.total_number_density(particle_2.pos.x, particle_2.pos.y, particle_2.pos.z); //We just need the lindhard screening length here, so the particular potential is not important let a: f64 = interactions::screening_length(Za, Zb, InteractionPotential::MOLIERE); let reduced_energy: f64 = LINDHARD_REDUCED_ENERGY_PREFACTOR*a*Mb/(Ma+Mb)/Za/Zb*E; let estimated_range_of_recoils = (reduced_energy.powf(0.3) + 0.1).powf(3.)/n/a/a; - if let Closest::SinglePoint(p2) = material.closest_point(particle_2.pos.x, particle_2.pos.y) { + //TODO: CHange this 2D + if let Closest::SinglePoint(p2) = material.closest_point(particle_2.pos.x, particle_2.pos.y, particle_2.pos.z) { let dx = p2.x() - particle_2.pos.x; let dy = p2.y() - particle_2.pos.y; let distance_to_surface = (dx*dx + dy*dy).sqrt(); @@ -200,7 +201,7 @@ pub fn determine_mfp_phi_impact_parameter(particle_1: &mut particle::Particle, m let y = particle_1.pos.y; let z = particle_1.pos.z; - let mut mfp = material.mfp(x, y); + let mut mfp = material.mfp(x, y, z); let mut phis_azimuthal = Vec::with_capacity(options.weak_collision_order + 1); let mut binary_collision_geometries = Vec::with_capacity(options.weak_collision_order + 1); @@ -214,11 +215,11 @@ pub fn determine_mfp_phi_impact_parameter(particle_1: &mut particle::Particle, m if options.high_energy_free_flight_paths { let Ma: f64 = particle_1.m; - let Mb: f64 = material.average_mass(x, y); + let Mb: f64 = material.average_mass(x, y, z); let Za: f64 = particle_1.Z; - let Zb: f64 = material.average_Z(x, y); - let n: &Vec = material.number_densities(x, y); - let ck: f64 = material.electronic_stopping_correction_factor(x, y); + let Zb: f64 = material.average_Z(x, y, z); + let n: &Vec = material.number_densities(x, y, z); + let ck: f64 = material.electronic_stopping_correction_factor(x, y, z); let E: f64 = particle_1.E; let Ec: f64 = particle_1.Ec; //We just need the Lindhard screening length here, so the particular potential is not important @@ -232,16 +233,16 @@ pub fn determine_mfp_phi_impact_parameter(particle_1: &mut particle::Particle, m //Free flight path formulation here described in SRIM textbook chapter 7, and Eckstein 1991 7.5.2 let ep = (reduced_energy*reduced_energy_min).sqrt(); let mut pmax = a/(ep + ep.sqrt() + 0.125*ep.powf(0.1)); - let mut ffp = 1./(material.total_number_density(x, y)*pmax*pmax*PI); + let mut ffp = 1./(material.total_number_density(x, y, z)*pmax*pmax*PI); let stopping_powers = material.electronic_stopping_cross_sections(particle_1, ElectronicStoppingMode::INTERPOLATED); - let delta_energy_electronic = stopping_powers.iter().zip(material.number_densities(x, y)).map(|(&se, number_density)| se*number_density).sum::()*ffp*ck; + let delta_energy_electronic = stopping_powers.iter().zip(material.number_densities(x, y, z)).map(|(&se, number_density)| se*number_density).sum::()*ffp*ck; //If losing too much energy, scale free-flight-path down //5 percent limit set in original TRIM paper, Biersack and Haggmark 1980 if delta_energy_electronic > 0.05*E { ffp *= 0.05*E/delta_energy_electronic; - pmax = (1./(material.total_number_density(x, y)*PI*ffp)).sqrt() + pmax = (1./(material.total_number_density(x, y, z)*PI*ffp)).sqrt() } //If free-flight-path less than the interatomic spacing, revert to solid model @@ -348,7 +349,7 @@ pub fn choose_collision_partner(particle_1: &particle::Particle, material: &mate let z_recoil: f64 = z + mfp*cosz + impact_parameter*(sinphi*cosy - cosphi*cosx*cosz)/sinx; //Choose recoil Z, M - let (species_index, Z_recoil, M_recoil, Ec_recoil, Es_recoil, interaction_index) = material.choose(x_recoil, y_recoil); + let (species_index, Z_recoil, M_recoil, Ec_recoil, Es_recoil, interaction_index) = material.choose(x_recoil, y_recoil, z_recoil); return (species_index, particle::Particle::new( @@ -417,12 +418,13 @@ pub fn update_particle_energy(particle_1: &mut particle::Particle, material: &ma let x = particle_1.pos.x; let y = particle_1.pos.y; + let z = particle_1.pos.z; - if material.inside(x, y) { + if material.inside(x, y, z) { let interaction_potential = options.interaction_potential[particle_1.interaction_index][material.interaction_index[strong_collision_index]]; let electronic_stopping_powers = material.electronic_stopping_cross_sections(particle_1, options.electronic_stopping_mode); - let n = material.number_densities(x, y); + let n = material.number_densities(x, y, z); let delta_energy = match options.electronic_stopping_mode { ElectronicStoppingMode::INTERPOLATED => electronic_stopping_powers.iter().zip(n).map(|(se, number_density)| se*number_density).collect::>().iter().sum::()*distance_traveled, diff --git a/src/main.rs b/src/main.rs index abd4aba..4ff9963 100644 --- a/src/main.rs +++ b/src/main.rs @@ -64,6 +64,7 @@ pub use crate::consts::*; pub use crate::structs::*; pub use crate::input::{Input, Options}; pub use crate::output::{OutputUnits}; +pub use crate::mesh::{Mesh, MeshElement}; fn main() { //Open and process input_file diff --git a/src/material.rs b/src/material.rs index 86ef326..176abd2 100644 --- a/src/material.rs +++ b/src/material.rs @@ -67,22 +67,22 @@ impl Material { } /// Gets concentrations of triangle that contains or is nearest to (x, y) - pub fn get_concentrations(&self, x: f64, y: f64) -> Vec { + pub fn get_concentrations(&self, x: f64, y: f64, z: f64) -> Vec { - let total_number_density: f64 = if self.inside(x, y) { - self.mesh_2d.get_total_density(x,y) + let total_number_density: f64 = if self.inside(x, y, z) { + self.mesh_2d.get_total_density(x, y, z) } else { - self.mesh_2d.nearest_cell_to(x, y).densities.iter().sum::() + self.mesh_2d.nearest_cell_to(x, y, z).densities.iter().sum::() }; - return self.mesh_2d.get_densities(x, y).iter().map(|&i| i / total_number_density).collect(); + return self.mesh_2d.get_densities(x, y, z).iter().map(|&i| i / total_number_density).collect(); } /// Gets cumulative concentrations of triangle that contains or is nearest to (x, y). /// Used for weighting the random choice of one of the constituent species. - pub fn get_cumulative_concentrations(&self, x: f64, y: f64) -> Vec { + pub fn get_cumulative_concentrations(&self, x: f64, y: f64, z: f64) -> Vec { let mut sum: f64 = 0.; - let concentrations = self.mesh_2d.get_concentrations(x, y); + let concentrations = self.mesh_2d.get_concentrations(x, y, z); let mut cumulative_concentrations = Vec::with_capacity(concentrations.len()); for concentration in concentrations { sum += concentration; @@ -92,52 +92,52 @@ impl Material { } /// Determines whether (x, y) is inside the material. - pub fn inside(&self, x: f64, y: f64) -> bool { - return self.mesh_2d.inside(x, y); + pub fn inside(&self, x: f64, y: f64, z: f64) -> bool { + return self.mesh_2d.inside(x, y, z); } /// Gets electronic stopping correction factor for LS and OR - pub fn electronic_stopping_correction_factor(&self, x: f64, y: f64) -> f64 { - if self.inside(x, y) { - return self.mesh_2d.get_ck(x, y); + pub fn electronic_stopping_correction_factor(&self, x: f64, y: f64, z: f64) -> f64 { + if self.inside(x, y, z) { + return self.mesh_2d.get_ck(x, y, z); } else { - return self.mesh_2d.nearest_cell_to(x, y).electronic_stopping_correction_factor; + return self.mesh_2d.nearest_cell_to(x, y, z).electronic_stopping_correction_factor; } } /// Determines the local mean free path from the formula sum(n(x, y))^(-1/3) - pub fn mfp(&self, x: f64, y: f64) -> f64 { - return self.total_number_density(x, y).powf(-1./3.); + pub fn mfp(&self, x: f64, y: f64, z: f64) -> f64 { + return self.total_number_density(x, y, z).powf(-1./3.); } /// Total number density of triangle that contains or is nearest to (x, y). - pub fn total_number_density(&self, x: f64, y: f64) -> f64 { - if self.inside(x, y) { - return self.mesh_2d.get_densities(x, y).iter().sum::(); + pub fn total_number_density(&self, x: f64, y: f64, z: f64) -> f64 { + if self.inside(x, y, z) { + return self.mesh_2d.get_densities(x, y, z).iter().sum::(); } else { - return self.mesh_2d.nearest_cell_to(x, y).densities.iter().sum::(); + return self.mesh_2d.nearest_cell_to(x, y, z).densities.iter().sum::(); } //return self.n.iter().sum::(); } /// Lists number density of each species of triangle that contains or is nearest to (x, y). - pub fn number_densities(&self, x: f64, y: f64) -> &Vec { - if self.inside(x, y) { - return &self.mesh_2d.get_densities(x, y); + pub fn number_densities(&self, x: f64, y: f64, z: f64) -> &Vec { + if self.inside(x, y, z) { + return &self.mesh_2d.get_densities(x, y, z); } else { - return &self.mesh_2d.nearest_cell_to(x, y).densities; + return &self.mesh_2d.nearest_cell_to(x, y, z).densities; } //return self.n.iter().sum::(); } /// Determines whether a point (x, y) is inside the energy barrier of the material. - pub fn inside_energy_barrier(&self, x: f64, y: f64) -> bool { + pub fn inside_energy_barrier(&self, x: f64, y: f64, z: f64) -> bool { - if self.mesh_2d.inside(x, y) { + if self.mesh_2d.inside(x, y, z) { true } else { - let nearest_cell = self.mesh_2d.nearest_cell_to(x, y); - let distance = nearest_cell.distance_to(x, y); + let nearest_cell = self.mesh_2d.nearest_cell_to(x, y, z); + let distance = nearest_cell.distance_to(x, y, z); distance < self.mesh_2d.energy_barrier_thickness } //let p = point!(x: x, y: y); @@ -145,42 +145,44 @@ impl Material { } /// Determines whether a point (x, y) is inside the simulation boundary. - pub fn inside_simulation_boundary(&self, x:f64, y: f64) -> bool { + pub fn inside_simulation_boundary(&self, x:f64, y: f64, z: f64) -> bool { let p = point!(x: x, y: y); return self.mesh_2d.simulation_boundary.contains(&p); } + //TODO: CHANGE THIS TO REMOVE IMPLIED 2D /// Finds the closest point on the material boundary to the point (x, y). - pub fn closest_point(&self, x: f64, y: f64) -> Closest { + pub fn closest_point(&self, x: f64, y: f64, z: f64) -> Closest { let p = point!(x: x, y: y); return self.mesh_2d.boundary.closest_point(&p); //return self.surface.closest_point(&p) } + //TODO: CHANGE THIS TO REMOVE IMPLIED 2D /// Finds the closest point on the simulation boundary to the point (x, y). - pub fn closest_point_on_simulation_surface(&self, x: f64, y: f64) -> Closest { + pub fn closest_point_on_simulation_surface(&self, x: f64, y: f64, z: f64) -> Closest { let p = point!(x: x, y: y); return self.mesh_2d.simulation_boundary.closest_point(&p) } /// Finds the average, concentration-weighted atomic number, Z_effective, of the triangle that contains or is nearest to (x, y). - pub fn average_Z(&self, x: f64, y: f64) -> f64 { - let concentrations = self.mesh_2d.get_concentrations(x, y); + pub fn average_Z(&self, x: f64, y: f64, z: f64) -> f64 { + let concentrations = self.mesh_2d.get_concentrations(x, y, z); return self.Z.iter().zip(concentrations).map(|(charge, concentration)| charge*concentration).collect::>().iter().sum(); //return self.Z.iter().sum::()/self.Z.len() as f64; } /// Finds the average, concentration-weighted atomic mass, m_effective, of the triangle that contains or is nearest to (x, y). - pub fn average_mass(&self, x: f64, y: f64) -> f64 { - let concentrations = self.mesh_2d.get_concentrations(x, y); + pub fn average_mass(&self, x: f64, y: f64, z: f64) -> f64 { + let concentrations = self.mesh_2d.get_concentrations(x, y, z); return self.m.iter().zip(concentrations).map(|(mass, concentration)| mass*concentration).collect::>().iter().sum(); //return self.m.iter().sum::()/self.m.len() as f64; } /// Finds the average, concentration-weighted bulk binding energy of the triangle that contains or is nearest to (x, y). - pub fn average_bulk_binding_energy(&self, x: f64, y: f64) -> f64 { + pub fn average_bulk_binding_energy(&self, x: f64, y: f64, z: f64) -> f64 { //returns average bulk binding energy - let concentrations = self.mesh_2d.get_concentrations(x, y); + let concentrations = self.mesh_2d.get_concentrations(x, y, z); return self.Eb.iter().zip(concentrations).map(|(bulk_binding_energy, concentration)| bulk_binding_energy*concentration).collect::>().iter().sum(); //return self.Eb.iter().sum::()/self.Eb.len() as f64; } @@ -190,8 +192,8 @@ impl Material { /// 1. INDIVIDUAL: the surface binding energies are set individually for each species, as Es. /// 2. TARGET: the surface binding energy is calculated as the local concentration-weighted average of the target surface binding energies, unless the particle has Es = 0, in which case it is 0. /// 3. AVERAGE: the surface binding energy is the average of the particle surface binding energy and the TARGET surface binding energy, unless either is 0 in which case it is 0. - pub fn actual_surface_binding_energy(&self, particle: &particle::Particle, x: f64, y: f64) -> f64 { - let concentrations = self.mesh_2d.get_concentrations(x, y); + pub fn actual_surface_binding_energy(&self, particle: &particle::Particle, x: f64, y: f64, z: f64) -> f64 { + let concentrations = self.mesh_2d.get_concentrations(x, y, z); match self.surface_binding_model { SurfaceBindingModel::INDIVIDUAL => particle.Es, @@ -227,9 +229,9 @@ impl Material { } ///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) -> (usize, f64, f64, f64, f64, usize) { + pub fn choose(&self, x: f64, y: f64, z: f64) -> (usize, f64, f64, f64, f64, usize) { let random_number: f64 = rand::random::(); - let cumulative_concentrations = self.get_cumulative_concentrations(x, y); + let cumulative_concentrations = self.get_cumulative_concentrations(x, y, z); for (component_index, cumulative_concentration) in cumulative_concentrations.iter().enumerate() { if random_number < *cumulative_concentration { @@ -253,9 +255,10 @@ impl Material { let x = particle_1.pos.x; let y = particle_1.pos.y; - let ck = self.electronic_stopping_correction_factor(x, y); + let z = particle_1.pos.y; + let ck = self.electronic_stopping_correction_factor(x, y, z); - for (n, Zb) in self.number_densities(x, y).iter().zip(&self.Z) { + for (n, Zb) in self.number_densities(x, y, z).iter().zip(&self.Z) { //let n = self.number_density(particle_1.pos.x, particle_1.pos.y); //let Zb = self.Z_eff(particle_1.pos.x, particle_1.pos.y); @@ -304,27 +307,30 @@ impl Material { pub fn surface_binding_energy(particle_1: &mut particle::Particle, material: &material::Material) { let x = particle_1.pos.x; let y = particle_1.pos.y; + let z = particle_1.pos.z; let x_old = particle_1.pos_old.x; let y_old = particle_1.pos_old.y; + let z_old = particle_1.pos_old.z; let cosx = particle_1.dir.x; let cosy = particle_1.dir.y; let cosz = particle_1.dir.z; let E = particle_1.E; //Actual surface binding energies - let Es = material.actual_surface_binding_energy(particle_1, x_old, y_old); + let Es = material.actual_surface_binding_energy(particle_1, x_old, y_old, z_old); let Ec = particle_1.Ec; - let inside_now = material.inside_energy_barrier(x, y); - let inside_old = material.inside_energy_barrier(x_old, y_old); + let inside_now = material.inside_energy_barrier(x, y, z); + let inside_old = material.inside_energy_barrier(x_old, y_old, z_old); let leaving = !inside_now & inside_old; let entering = inside_now & !inside_old; + //TODO: Remove explicit 2D dependence if entering { if particle_1.backreflected { particle_1.backreflected = false; - } else if let Closest::SinglePoint(p2) = material.closest_point(x, y) { + } else if let Closest::SinglePoint(p2) = material.closest_point(x, y, z) { let dx = p2.x() - x; let dy = p2.y() - y; let mag = (dx*dx + dy*dy).sqrt(); @@ -343,7 +349,7 @@ pub fn surface_binding_energy(particle_1: &mut particle::Particle, material: &ma } if leaving { - if let Closest::SinglePoint(p2) = material.closest_point(x, y) { + if let Closest::SinglePoint(p2) = material.closest_point(x, y, z) { let dx = p2.x() - x; let dy = p2.y() - y; @@ -384,15 +390,16 @@ pub fn surface_binding_energy(particle_1: &mut particle::Particle, material: &ma pub fn boundary_condition_2D_planar(particle_1: &mut particle::Particle, material: &material::Material) { let x = particle_1.pos.x; let y = particle_1.pos.y; + let z = particle_1.pos.z; let E = particle_1.E; let Ec = particle_1.Ec; - if !material.inside_simulation_boundary(x, y) { + if !material.inside_simulation_boundary(x, y, z) { particle_1.left = true; particle_1.add_trajectory(); } - if (E < Ec) & !particle_1.left & material.inside_energy_barrier(x, y) { + if (E < Ec) & !particle_1.left & material.inside_energy_barrier(x, y, z) { particle_1.stopped = true; particle_1.add_trajectory(); } diff --git a/src/mesh.rs b/src/mesh.rs index af1ded2..23b83fb 100644 --- a/src/mesh.rs +++ b/src/mesh.rs @@ -22,6 +22,21 @@ pub struct Mesh2D { pub simulation_boundary: Polygon, pub energy_barrier_thickness: f64 } + +pub trait Mesh { + fn get_densities(&self, x: f64, y: f64, z: f64) -> &Vec; + fn get_ck(&self, x: f64, y: f64, z: f64) -> f64; + fn get_total_density(&self, x: f64, y: f64, z: f64) -> f64; + fn get_concentrations(&self, x: f64, y: f64, z: f64) -> &Vec; + fn nearest_cell_to(&self, x: f64, y: f64, z: f64) -> &Cell2D; + fn inside(&self, x: f64, y: f64, z: f64) -> bool; +} + +pub trait MeshElement { + fn contains(&self, x: f64, y: f64, z: f64) -> bool; + fn distance_to(&self, x: f64, y: f64, z: f64) -> f64; +} + impl Mesh2D { /// Constructor for Mesh2D object from mesh_2d_input. pub fn new(mesh_2d_input: Mesh2DInput) -> Mesh2D { @@ -98,11 +113,12 @@ impl Mesh2D { energy_barrier_thickness: energy_barrier_thickness, } } - +} +impl Mesh for Mesh2D { /// Find the number densities of the triangle that contains or is nearest to (x, y). - pub fn get_densities(&self, x: f64, y: f64) -> &Vec { + fn get_densities(&self, x: f64, y: f64, z: f64) -> &Vec { for cell in &self.mesh { - if cell.contains(x, y) { + if cell.contains(x, y, z) { return &cell.densities; } } @@ -110,9 +126,9 @@ impl Mesh2D { } /// Find the number densities of the triangle that contains or is nearest to (x, y). - pub fn get_ck(&self, x: f64, y: f64) -> f64 { + fn get_ck(&self, x: f64, y: f64, z: f64) -> f64 { for cell in &self.mesh { - if cell.contains(x, y) { + if cell.contains(x, y, z) { return cell.electronic_stopping_correction_factor; } } @@ -120,9 +136,9 @@ impl Mesh2D { } /// Determine the total number density of the triangle that contains or is nearest to (x, y). - pub fn get_total_density(&self, x: f64, y: f64) -> f64 { + fn get_total_density(&self, x: f64, y: f64, z: f64) -> f64 { for cell in &self.mesh { - if cell.contains(x, y) { + if cell.contains(x, y, z) { return cell.densities.iter().sum::(); } } @@ -130,32 +146,32 @@ impl Mesh2D { } /// Find the concentrations of the triangle that contains or is nearest to (x, y). - pub fn get_concentrations(&self, x: f64, y: f64) -> &Vec { - if self.inside(x, y) { + fn get_concentrations(&self, x: f64, y: f64, z: f64) -> &Vec { + if self.inside(x, y, z) { for cell in &self.mesh { - if cell.contains(x, y) { + if cell.contains(x, y, z) { return &cell.concentrations; } } panic!("Geometry error: method inside() is returning true for points outside all cells. Check boundary points.") } else { - return &self.nearest_cell_to(x, y).concentrations; + return &self.nearest_cell_to(x, y, z).concentrations; } } /// Determines whether the point (x, y) is inside the mesh. - pub fn inside(&self, x: f64, y: f64) -> bool { + fn inside(&self, x: f64, y: f64, z: f64) -> bool { self.boundary.contains(&Point::new(x, y)) } /// Finds the cell that is nearest to (x, y). - pub fn nearest_cell_to(&self, x: f64, y: f64) -> &Cell2D { + fn nearest_cell_to(&self, x: f64, y: f64, z: f64) -> &Cell2D { let mut min_distance: f64 = std::f64::MAX; let mut index: usize = 0; for (cell_index, cell) in self.mesh.iter().enumerate() { - let distance_to = cell.distance_to(x, y); + let distance_to = cell.distance_to(x, y, z); if distance_to < min_distance { min_distance = distance_to; index = cell_index; @@ -166,6 +182,7 @@ impl Mesh2D { } } + /// A mesh cell that contains a triangle and the local number densities and concentrations. pub struct Cell2D { triangle: Triangle2D, @@ -183,14 +200,16 @@ impl Cell2D { electronic_stopping_correction_factor } } +} +impl MeshElement for Cell2D { /// Determines whether this cell contains the point (x, y). - pub fn contains(&self, x: f64, y: f64) -> bool { + fn contains(&self, x: f64, y: f64, z: f64) -> bool { self.triangle.contains(x, y) } /// Computes the shortest distance between this cell and the point (x, y). - pub fn distance_to(&self, x: f64, y: f64) -> f64 { + fn distance_to(&self, x: f64, y: f64, z: f64) -> f64 { self.triangle.distance_to(x, y) } } diff --git a/src/tests.rs b/src/tests.rs index 53ce3b6..dbcbcde 100644 --- a/src/tests.rs +++ b/src/tests.rs @@ -50,35 +50,35 @@ fn test_surface_binding_energy_barrier() { particle_1.pos_old.x = -500.*ANGSTROM; particle_1.pos_old.y = 0.; - let inside = material_1.inside(particle_1.pos.x, particle_1.pos.y); - let inside_old = material_1.inside(particle_1.pos_old.x, particle_1.pos_old.y); + let inside = material_1.inside(particle_1.pos.x, particle_1.pos.y, particle_1.pos.z); + let inside_old = material_1.inside(particle_1.pos_old.x, particle_1.pos_old.y, particle_1.pos.z); //println!("{} {}", inside, inside_old); assert!(inside); assert!(!inside_old); //Test concentration-dependent surface binding energy - let surface_binding_energy = material_1.actual_surface_binding_energy(&particle_1, particle_1.pos.x, particle_1.pos.y); + let surface_binding_energy = material_1.actual_surface_binding_energy(&particle_1, particle_1.pos.x, particle_1.pos.y, particle_1.pos.z); assert!(approx_eq!(f64, surface_binding_energy/EV, (2. + 4.)/2., epsilon=1E-12)); //println!("sbv: {}", surface_binding_energy/EV); //Test leftmost boundary - assert!(material_1.inside_energy_barrier(500.*ANGSTROM, 0.)); - assert!(material_1.inside_energy_barrier(-5.*ANGSTROM, 0.)); - assert!(!material_1.inside_energy_barrier(-15.*ANGSTROM, 0.)); + assert!(material_1.inside_energy_barrier(500.*ANGSTROM, 0., 0.)); + assert!(material_1.inside_energy_barrier(-5.*ANGSTROM, 0., 0.)); + assert!(!material_1.inside_energy_barrier(-15.*ANGSTROM, 0., 0.)); //Test top boundary - assert!(material_1.inside_energy_barrier(500.*ANGSTROM, 0.)); - assert!(material_1.inside_energy_barrier(500.*ANGSTROM, 505.*ANGSTROM)); - assert!(!material_1.inside_energy_barrier(500.*ANGSTROM, 515.*ANGSTROM)); + assert!(material_1.inside_energy_barrier(500.*ANGSTROM, 0., 0.)); + assert!(material_1.inside_energy_barrier(500.*ANGSTROM, 505.*ANGSTROM, 0.)); + assert!(!material_1.inside_energy_barrier(500.*ANGSTROM, 515.*ANGSTROM, 0.)); //Test bottom boundary - assert!(material_1.inside_energy_barrier(500.*ANGSTROM, -505.*ANGSTROM)); - assert!(!material_1.inside_energy_barrier(500.*ANGSTROM, -515.*ANGSTROM)); + assert!(material_1.inside_energy_barrier(500.*ANGSTROM, -505.*ANGSTROM, 0.)); + assert!(!material_1.inside_energy_barrier(500.*ANGSTROM, -515.*ANGSTROM, 0.)); //Test rightmost boundary - assert!(material_1.inside_energy_barrier(1005.*ANGSTROM, 0.)); - assert!(!material_1.inside_energy_barrier(1015.*ANGSTROM, 0.)); + assert!(material_1.inside_energy_barrier(1005.*ANGSTROM, 0., 0.)); + assert!(!material_1.inside_energy_barrier(1015.*ANGSTROM, 0., 0.)); } #[test] @@ -281,7 +281,7 @@ fn test_momentum_conservation() { println!("Initial Energies: {} eV {} eV", particle_1.E/EV, particle_2.E/EV); //Energy transfer to recoil - particle_2.E = binary_collision_result.recoil_energy - material_1.average_bulk_binding_energy(particle_2.pos.x, particle_2.pos.y); + particle_2.E = binary_collision_result.recoil_energy - material_1.average_bulk_binding_energy(particle_2.pos.x, particle_2.pos.y, particle_2.pos.z); //Rotate particle 1, 2 by lab frame scattering angles particle::rotate_particle(&mut particle_1, binary_collision_result.psi, From 4acd56bd326dedf537632bfa40ec69d47aadca8a Mon Sep 17 00:00:00 2001 From: Jon Drobny Date: Sat, 27 Feb 2021 12:54:49 -0600 Subject: [PATCH 02/14] This version has switched all methods with Material (with the exception of the default input) to Material or Material where necessary. From now on, any struct that implements the Mesh trait can be used to create a Material. Neat! --- src/bca.rs | 8 ++++---- src/input.rs | 4 ++-- src/material.rs | 22 +++++++++++----------- src/mesh.rs | 15 +++++++++++++++ src/tests.rs | 4 ++-- 5 files changed, 34 insertions(+), 19 deletions(-) diff --git a/src/bca.rs b/src/bca.rs index 624e6dc..8a152e0 100644 --- a/src/bca.rs +++ b/src/bca.rs @@ -63,7 +63,7 @@ impl fmt::Display for BinaryCollisionResult { } /// This function takes a single particle, a material, and an options object and runs a binary-collision-approximation trajectory for that particle in that material, producing a final particle list that consists of the original ion and any material particles displaced thereby. -pub fn single_ion_bca(particle: particle::Particle, material: &material::Material, options: &Options) -> Vec { +pub fn single_ion_bca(particle: particle::Particle, material: &material::Material, options: &Options) -> Vec { let mut particles: Vec = Vec::new(); particles.push(particle); @@ -195,7 +195,7 @@ pub fn single_ion_bca(particle: particle::Particle, material: &material::Materia /// For a particle in a material, determine the mean free path and choose the azimuthal angle and impact parameter. /// The mean free path can be exponentially distributed (gaseous) or constant (amorphous solid/liquid). Azimuthal angles are chosen uniformly. Impact parameters are chosen for collision partners distributed uniformly on a disk of density-dependent radius. -pub fn determine_mfp_phi_impact_parameter(particle_1: &mut particle::Particle, material: &material::Material, options: &Options) -> Vec { +pub fn determine_mfp_phi_impact_parameter(particle_1: &mut particle::Particle, material: &material::Material, options: &Options) -> Vec { let x = particle_1.pos.x; let y = particle_1.pos.y; @@ -326,7 +326,7 @@ pub fn determine_mfp_phi_impact_parameter(particle_1: &mut particle::Particle, m } /// For a particle in a material, and for a particular binary collision geometry, choose a species for the collision partner. -pub fn choose_collision_partner(particle_1: &particle::Particle, material: &material::Material, binary_collision_geometry: &BinaryCollisionGeometry, options: &Options) -> (usize, particle::Particle) { +pub fn choose_collision_partner(particle_1: &particle::Particle, material: &material::Material, binary_collision_geometry: &BinaryCollisionGeometry, options: &Options) -> (usize, particle::Particle) { let x = particle_1.pos.x; let y = particle_1.pos.y; let z = particle_1.pos.z; @@ -406,7 +406,7 @@ fn distance_of_closest_approach(particle_1: &particle::Particle, particle_2: &pa } /// Update a particle's energy following nuclear and electronic interactions of a single BCA step. -pub fn update_particle_energy(particle_1: &mut particle::Particle, material: &material::Material, distance_traveled: f64, +pub fn update_particle_energy(particle_1: &mut particle::Particle, material: &material::Material, distance_traveled: f64, recoil_energy: f64, x0: f64, strong_collision_Z: f64, strong_collision_index: usize, options: &Options) { //If particle energy drops below zero before electronic stopping calcualtion, it produces NaNs diff --git a/src/input.rs b/src/input.rs index de4d91a..8b1db8b 100644 --- a/src/input.rs +++ b/src/input.rs @@ -72,7 +72,7 @@ pub struct Options { pub z_num: usize, } -pub fn input() -> (Vec, material::Material, Options, OutputUnits){ +pub fn input() -> (Vec, material::Material, Options, OutputUnits){ let args: Vec = env::args().collect(); @@ -95,7 +95,7 @@ pub fn input() -> (Vec, material::Material, Options, Ou .expect("Input error: failed to parse TOML file. Check that all floats have terminal 0s and that there are no missing/extra delimiters."); //Unpack toml information into structs - let material = material::Material::new(input.material_parameters, input.mesh_2d_input); + let material: material::Material = material::Material::::new(input.material_parameters, input.mesh_2d_input); let options = input.options; let particle_parameters = input.particle_parameters; diff --git a/src/material.rs b/src/material.rs index 176abd2..9d9cbe5 100644 --- a/src/material.rs +++ b/src/material.rs @@ -19,20 +19,20 @@ pub struct MaterialParameters { } /// Material in rustbca. Includes the material properties and the mesh that defines the material geometry. -pub struct Material { +pub struct Material { pub m: Vec, pub Z: Vec, pub Eb: Vec, pub Es: Vec, pub Ec: Vec, pub interaction_index: Vec, - pub mesh_2d: mesh::Mesh2D, + pub mesh_2d: Box, pub surface_binding_model: SurfaceBindingModel } -impl Material { +impl Material { /// Constructs a new material object from a material parameters object and a mesh_2d_input object. - pub fn new(material_parameters: MaterialParameters, mesh_2d_input: mesh::Mesh2DInput) -> Material { + pub fn new(material_parameters: MaterialParameters, mesh_2d_input: mesh::Mesh2DInput) -> Material { let energy_unit: f64 = match material_parameters.energy_unit.as_str() { "EV" => EV, @@ -61,7 +61,7 @@ impl Material { Es: material_parameters.Es.iter().map(|&i| i*energy_unit).collect(), Ec: material_parameters.Ec.iter().map(|&i| i*energy_unit).collect(), interaction_index: material_parameters.interaction_index, - mesh_2d: mesh::Mesh2D::new(mesh_2d_input), + mesh_2d: Box::new(mesh::Mesh2D::new(mesh_2d_input)), surface_binding_model: material_parameters.surface_binding_model, } } @@ -138,7 +138,7 @@ impl Material { } else { let nearest_cell = self.mesh_2d.nearest_cell_to(x, y, z); let distance = nearest_cell.distance_to(x, y, z); - distance < self.mesh_2d.energy_barrier_thickness + distance < self.mesh_2d.get_energy_barrier_thickness() } //let p = point!(x: x, y: y); //return self.energy_surface.contains(&p); @@ -147,14 +147,14 @@ impl Material { /// Determines whether a point (x, y) is inside the simulation boundary. pub fn inside_simulation_boundary(&self, x:f64, y: f64, z: f64) -> bool { let p = point!(x: x, y: y); - return self.mesh_2d.simulation_boundary.contains(&p); + return self.mesh_2d.get_simulation_boundary().contains(&p); } //TODO: CHANGE THIS TO REMOVE IMPLIED 2D /// Finds the closest point on the material boundary to the point (x, y). pub fn closest_point(&self, x: f64, y: f64, z: f64) -> Closest { let p = point!(x: x, y: y); - return self.mesh_2d.boundary.closest_point(&p); + return self.mesh_2d.get_boundary().closest_point(&p); //return self.surface.closest_point(&p) } @@ -162,7 +162,7 @@ impl Material { /// Finds the closest point on the simulation boundary to the point (x, y). pub fn closest_point_on_simulation_surface(&self, x: f64, y: f64, z: f64) -> Closest { let p = point!(x: x, y: y); - return self.mesh_2d.simulation_boundary.closest_point(&p) + return self.mesh_2d.get_simulation_boundary().closest_point(&p) } /// Finds the average, concentration-weighted atomic number, Z_effective, of the triangle that contains or is nearest to (x, y). @@ -304,7 +304,7 @@ impl Material { /// Calculate the effects of the planar surface binding potential of a material on a particle. /// These effects include surface reflection and refraction of particles with non-zero surface binding energies. -pub fn surface_binding_energy(particle_1: &mut particle::Particle, material: &material::Material) { +pub fn surface_binding_energy(particle_1: &mut particle::Particle, material: &material::Material) { let x = particle_1.pos.x; let y = particle_1.pos.y; let z = particle_1.pos.z; @@ -387,7 +387,7 @@ pub fn surface_binding_energy(particle_1: &mut particle::Particle, material: &ma } /// Apply the boundary conditions of a material on a particle, including stopping, leaving, and reflection/refraction by/through the surface binding potential. -pub fn boundary_condition_2D_planar(particle_1: &mut particle::Particle, material: &material::Material) { +pub fn boundary_condition_2D_planar(particle_1: &mut particle::Particle, material: &material::Material) { let x = particle_1.pos.x; let y = particle_1.pos.y; let z = particle_1.pos.z; diff --git a/src/mesh.rs b/src/mesh.rs index 23b83fb..0fee416 100644 --- a/src/mesh.rs +++ b/src/mesh.rs @@ -30,6 +30,9 @@ pub trait Mesh { fn get_concentrations(&self, x: f64, y: f64, z: f64) -> &Vec; fn nearest_cell_to(&self, x: f64, y: f64, z: f64) -> &Cell2D; fn inside(&self, x: f64, y: f64, z: f64) -> bool; + fn get_energy_barrier_thickness(&self) -> f64; + fn get_simulation_boundary(&self) -> &Polygon; + fn get_boundary(&self) -> &Polygon; } pub trait MeshElement { @@ -115,6 +118,18 @@ impl Mesh2D { } } impl Mesh for Mesh2D { + fn get_energy_barrier_thickness(&self) -> f64 { + self.energy_barrier_thickness + } + + fn get_simulation_boundary(&self) -> &Polygon { + &self.simulation_boundary + } + + fn get_boundary(&self) -> &Polygon { + &self.boundary + } + /// Find the number densities of the triangle that contains or is nearest to (x, y). fn get_densities(&self, x: f64, y: f64, z: f64) -> &Vec { for cell in &self.mesh { diff --git a/src/tests.rs b/src/tests.rs index dbcbcde..5db578b 100644 --- a/src/tests.rs +++ b/src/tests.rs @@ -42,7 +42,7 @@ fn test_surface_binding_energy_barrier() { electronic_stopping_correction_factors: vec![0.0, 0.0], }; - let material_1 = material::Material::new(material_parameters, mesh_2d_input); + let material_1: material::Material = material::Material::::new(material_parameters, mesh_2d_input); particle_1.pos.x = 500.*ANGSTROM; particle_1.pos.y = 0.; @@ -225,7 +225,7 @@ fn test_momentum_conservation() { energy_barrier_thickness: 0., }; - let material_1 = material::Material::new(material_parameters, mesh_2d_input); + let material_1: material::Material = material::Material::::new(material_parameters, mesh_2d_input); for high_energy_free_flight_paths in vec![true, false] { for potential in vec![InteractionPotential::KR_C, InteractionPotential::MOLIERE, InteractionPotential::ZBL, InteractionPotential::LENZ_JENSEN] { From 445a101bf3a14697c16a6908a25076f7c9cc43e7 Mon Sep 17 00:00:00 2001 From: Jon Drobny Date: Sat, 27 Feb 2021 12:58:34 -0600 Subject: [PATCH 03/14] Nomenclature change - in Material, mesh_2d has become geometry. --- src/material.rs | 46 +++++++++++++++++++++++----------------------- 1 file changed, 23 insertions(+), 23 deletions(-) diff --git a/src/material.rs b/src/material.rs index 9d9cbe5..e8671c0 100644 --- a/src/material.rs +++ b/src/material.rs @@ -26,7 +26,7 @@ pub struct Material { pub Es: Vec, pub Ec: Vec, pub interaction_index: Vec, - pub mesh_2d: Box, + pub geometry: Box, pub surface_binding_model: SurfaceBindingModel } @@ -61,7 +61,7 @@ impl Material { Es: material_parameters.Es.iter().map(|&i| i*energy_unit).collect(), Ec: material_parameters.Ec.iter().map(|&i| i*energy_unit).collect(), interaction_index: material_parameters.interaction_index, - mesh_2d: Box::new(mesh::Mesh2D::new(mesh_2d_input)), + geometry: Box::new(mesh::Mesh2D::new(mesh_2d_input)), surface_binding_model: material_parameters.surface_binding_model, } } @@ -70,19 +70,19 @@ impl Material { pub fn get_concentrations(&self, x: f64, y: f64, z: f64) -> Vec { let total_number_density: f64 = if self.inside(x, y, z) { - self.mesh_2d.get_total_density(x, y, z) + self.geometry.get_total_density(x, y, z) } else { - self.mesh_2d.nearest_cell_to(x, y, z).densities.iter().sum::() + self.geometry.nearest_cell_to(x, y, z).densities.iter().sum::() }; - return self.mesh_2d.get_densities(x, y, z).iter().map(|&i| i / total_number_density).collect(); + return self.geometry.get_densities(x, y, z).iter().map(|&i| i / total_number_density).collect(); } /// Gets cumulative concentrations of triangle that contains or is nearest to (x, y). /// Used for weighting the random choice of one of the constituent species. pub fn get_cumulative_concentrations(&self, x: f64, y: f64, z: f64) -> Vec { let mut sum: f64 = 0.; - let concentrations = self.mesh_2d.get_concentrations(x, y, z); + let concentrations = self.geometry.get_concentrations(x, y, z); let mut cumulative_concentrations = Vec::with_capacity(concentrations.len()); for concentration in concentrations { sum += concentration; @@ -93,15 +93,15 @@ impl Material { /// Determines whether (x, y) is inside the material. pub fn inside(&self, x: f64, y: f64, z: f64) -> bool { - return self.mesh_2d.inside(x, y, z); + return self.geometry.inside(x, y, z); } /// Gets electronic stopping correction factor for LS and OR pub fn electronic_stopping_correction_factor(&self, x: f64, y: f64, z: f64) -> f64 { if self.inside(x, y, z) { - return self.mesh_2d.get_ck(x, y, z); + return self.geometry.get_ck(x, y, z); } else { - return self.mesh_2d.nearest_cell_to(x, y, z).electronic_stopping_correction_factor; + return self.geometry.nearest_cell_to(x, y, z).electronic_stopping_correction_factor; } } @@ -113,9 +113,9 @@ impl Material { /// Total number density of triangle that contains or is nearest to (x, y). pub fn total_number_density(&self, x: f64, y: f64, z: f64) -> f64 { if self.inside(x, y, z) { - return self.mesh_2d.get_densities(x, y, z).iter().sum::(); + return self.geometry.get_densities(x, y, z).iter().sum::(); } else { - return self.mesh_2d.nearest_cell_to(x, y, z).densities.iter().sum::(); + return self.geometry.nearest_cell_to(x, y, z).densities.iter().sum::(); } //return self.n.iter().sum::(); } @@ -123,9 +123,9 @@ impl Material { /// Lists number density of each species of triangle that contains or is nearest to (x, y). pub fn number_densities(&self, x: f64, y: f64, z: f64) -> &Vec { if self.inside(x, y, z) { - return &self.mesh_2d.get_densities(x, y, z); + return &self.geometry.get_densities(x, y, z); } else { - return &self.mesh_2d.nearest_cell_to(x, y, z).densities; + return &self.geometry.nearest_cell_to(x, y, z).densities; } //return self.n.iter().sum::(); } @@ -133,12 +133,12 @@ impl Material { /// Determines whether a point (x, y) is inside the energy barrier of the material. pub fn inside_energy_barrier(&self, x: f64, y: f64, z: f64) -> bool { - if self.mesh_2d.inside(x, y, z) { + if self.geometry.inside(x, y, z) { true } else { - let nearest_cell = self.mesh_2d.nearest_cell_to(x, y, z); + let nearest_cell = self.geometry.nearest_cell_to(x, y, z); let distance = nearest_cell.distance_to(x, y, z); - distance < self.mesh_2d.get_energy_barrier_thickness() + distance < self.geometry.get_energy_barrier_thickness() } //let p = point!(x: x, y: y); //return self.energy_surface.contains(&p); @@ -147,14 +147,14 @@ impl Material { /// Determines whether a point (x, y) is inside the simulation boundary. pub fn inside_simulation_boundary(&self, x:f64, y: f64, z: f64) -> bool { let p = point!(x: x, y: y); - return self.mesh_2d.get_simulation_boundary().contains(&p); + return self.geometry.get_simulation_boundary().contains(&p); } //TODO: CHANGE THIS TO REMOVE IMPLIED 2D /// Finds the closest point on the material boundary to the point (x, y). pub fn closest_point(&self, x: f64, y: f64, z: f64) -> Closest { let p = point!(x: x, y: y); - return self.mesh_2d.get_boundary().closest_point(&p); + return self.geometry.get_boundary().closest_point(&p); //return self.surface.closest_point(&p) } @@ -162,19 +162,19 @@ impl Material { /// Finds the closest point on the simulation boundary to the point (x, y). pub fn closest_point_on_simulation_surface(&self, x: f64, y: f64, z: f64) -> Closest { let p = point!(x: x, y: y); - return self.mesh_2d.get_simulation_boundary().closest_point(&p) + return self.geometry.get_simulation_boundary().closest_point(&p) } /// Finds the average, concentration-weighted atomic number, Z_effective, of the triangle that contains or is nearest to (x, y). pub fn average_Z(&self, x: f64, y: f64, z: f64) -> f64 { - let concentrations = self.mesh_2d.get_concentrations(x, y, z); + let concentrations = self.geometry.get_concentrations(x, y, z); return self.Z.iter().zip(concentrations).map(|(charge, concentration)| charge*concentration).collect::>().iter().sum(); //return self.Z.iter().sum::()/self.Z.len() as f64; } /// Finds the average, concentration-weighted atomic mass, m_effective, of the triangle that contains or is nearest to (x, y). pub fn average_mass(&self, x: f64, y: f64, z: f64) -> f64 { - let concentrations = self.mesh_2d.get_concentrations(x, y, z); + let concentrations = self.geometry.get_concentrations(x, y, z); return self.m.iter().zip(concentrations).map(|(mass, concentration)| mass*concentration).collect::>().iter().sum(); //return self.m.iter().sum::()/self.m.len() as f64; } @@ -182,7 +182,7 @@ impl Material { /// Finds the average, concentration-weighted bulk binding energy of the triangle that contains or is nearest to (x, y). pub fn average_bulk_binding_energy(&self, x: f64, y: f64, z: f64) -> f64 { //returns average bulk binding energy - let concentrations = self.mesh_2d.get_concentrations(x, y, z); + let concentrations = self.geometry.get_concentrations(x, y, z); return self.Eb.iter().zip(concentrations).map(|(bulk_binding_energy, concentration)| bulk_binding_energy*concentration).collect::>().iter().sum(); //return self.Eb.iter().sum::()/self.Eb.len() as f64; } @@ -193,7 +193,7 @@ impl Material { /// 2. TARGET: the surface binding energy is calculated as the local concentration-weighted average of the target surface binding energies, unless the particle has Es = 0, in which case it is 0. /// 3. AVERAGE: the surface binding energy is the average of the particle surface binding energy and the TARGET surface binding energy, unless either is 0 in which case it is 0. pub fn actual_surface_binding_energy(&self, particle: &particle::Particle, x: f64, y: f64, z: f64) -> f64 { - let concentrations = self.mesh_2d.get_concentrations(x, y, z); + let concentrations = self.geometry.get_concentrations(x, y, z); match self.surface_binding_model { SurfaceBindingModel::INDIVIDUAL => particle.Es, From ac0ea9d0448f3559a3dbb51d26c544fc0e77e487 Mon Sep 17 00:00:00 2001 From: Jon Drobny Date: Sat, 27 Feb 2021 13:06:54 -0600 Subject: [PATCH 04/14] Replaced Mesh nomenclature with the more accurate Geometry. --- src/bca.rs | 8 ++++---- src/main.rs | 2 +- src/material.rs | 9 +++++---- src/mesh.rs | 24 ++++++++++++++++++++---- 4 files changed, 30 insertions(+), 13 deletions(-) diff --git a/src/bca.rs b/src/bca.rs index 8a152e0..ee3cd20 100644 --- a/src/bca.rs +++ b/src/bca.rs @@ -63,7 +63,7 @@ impl fmt::Display for BinaryCollisionResult { } /// This function takes a single particle, a material, and an options object and runs a binary-collision-approximation trajectory for that particle in that material, producing a final particle list that consists of the original ion and any material particles displaced thereby. -pub fn single_ion_bca(particle: particle::Particle, material: &material::Material, options: &Options) -> Vec { +pub fn single_ion_bca(particle: particle::Particle, material: &material::Material, options: &Options) -> Vec { let mut particles: Vec = Vec::new(); particles.push(particle); @@ -195,7 +195,7 @@ pub fn single_ion_bca(particle: particle::Particle, material: &material /// For a particle in a material, determine the mean free path and choose the azimuthal angle and impact parameter. /// The mean free path can be exponentially distributed (gaseous) or constant (amorphous solid/liquid). Azimuthal angles are chosen uniformly. Impact parameters are chosen for collision partners distributed uniformly on a disk of density-dependent radius. -pub fn determine_mfp_phi_impact_parameter(particle_1: &mut particle::Particle, material: &material::Material, options: &Options) -> Vec { +pub fn determine_mfp_phi_impact_parameter(particle_1: &mut particle::Particle, material: &material::Material, options: &Options) -> Vec { let x = particle_1.pos.x; let y = particle_1.pos.y; @@ -326,7 +326,7 @@ pub fn determine_mfp_phi_impact_parameter(particle_1: &mut particle::Pa } /// For a particle in a material, and for a particular binary collision geometry, choose a species for the collision partner. -pub fn choose_collision_partner(particle_1: &particle::Particle, material: &material::Material, binary_collision_geometry: &BinaryCollisionGeometry, options: &Options) -> (usize, particle::Particle) { +pub fn choose_collision_partner(particle_1: &particle::Particle, material: &material::Material, binary_collision_geometry: &BinaryCollisionGeometry, options: &Options) -> (usize, particle::Particle) { let x = particle_1.pos.x; let y = particle_1.pos.y; let z = particle_1.pos.z; @@ -406,7 +406,7 @@ fn distance_of_closest_approach(particle_1: &particle::Particle, particle_2: &pa } /// Update a particle's energy following nuclear and electronic interactions of a single BCA step. -pub fn update_particle_energy(particle_1: &mut particle::Particle, material: &material::Material, distance_traveled: f64, +pub fn update_particle_energy(particle_1: &mut particle::Particle, material: &material::Material, distance_traveled: f64, recoil_energy: f64, x0: f64, strong_collision_Z: f64, strong_collision_index: usize, options: &Options) { //If particle energy drops below zero before electronic stopping calcualtion, it produces NaNs diff --git a/src/main.rs b/src/main.rs index 4ff9963..a1c88d4 100644 --- a/src/main.rs +++ b/src/main.rs @@ -64,7 +64,7 @@ pub use crate::consts::*; pub use crate::structs::*; pub use crate::input::{Input, Options}; pub use crate::output::{OutputUnits}; -pub use crate::mesh::{Mesh, MeshElement}; +pub use crate::mesh::{Geometry, GeometryElement}; fn main() { //Open and process input_file diff --git a/src/material.rs b/src/material.rs index e8671c0..c9c0e85 100644 --- a/src/material.rs +++ b/src/material.rs @@ -19,7 +19,7 @@ pub struct MaterialParameters { } /// Material in rustbca. Includes the material properties and the mesh that defines the material geometry. -pub struct Material { +pub struct Material { pub m: Vec, pub Z: Vec, pub Eb: Vec, @@ -30,7 +30,8 @@ pub struct Material { pub surface_binding_model: SurfaceBindingModel } -impl Material { +impl Material { + /// Constructs a new material object from a material parameters object and a mesh_2d_input object. pub fn new(material_parameters: MaterialParameters, mesh_2d_input: mesh::Mesh2DInput) -> Material { @@ -304,7 +305,7 @@ impl Material { /// Calculate the effects of the planar surface binding potential of a material on a particle. /// These effects include surface reflection and refraction of particles with non-zero surface binding energies. -pub fn surface_binding_energy(particle_1: &mut particle::Particle, material: &material::Material) { +pub fn surface_binding_energy(particle_1: &mut particle::Particle, material: &material::Material) { let x = particle_1.pos.x; let y = particle_1.pos.y; let z = particle_1.pos.z; @@ -387,7 +388,7 @@ pub fn surface_binding_energy(particle_1: &mut particle::Particle, mate } /// Apply the boundary conditions of a material on a particle, including stopping, leaving, and reflection/refraction by/through the surface binding potential. -pub fn boundary_condition_2D_planar(particle_1: &mut particle::Particle, material: &material::Material) { +pub fn boundary_condition_2D_planar(particle_1: &mut particle::Particle, material: &material::Material) { let x = particle_1.pos.x; let y = particle_1.pos.y; let z = particle_1.pos.z; diff --git a/src/mesh.rs b/src/mesh.rs index 0fee416..6bb338e 100644 --- a/src/mesh.rs +++ b/src/mesh.rs @@ -23,7 +23,8 @@ pub struct Mesh2D { pub energy_barrier_thickness: f64 } -pub trait Mesh { +///Trait for a Geometry object +pub trait Geometry { fn get_densities(&self, x: f64, y: f64, z: f64) -> &Vec; fn get_ck(&self, x: f64, y: f64, z: f64) -> f64; fn get_total_density(&self, x: f64, y: f64, z: f64) -> f64; @@ -35,9 +36,12 @@ pub trait Mesh { fn get_boundary(&self) -> &Polygon; } -pub trait MeshElement { +pub trait GeometryElement { fn contains(&self, x: f64, y: f64, z: f64) -> bool; fn distance_to(&self, x: f64, y: f64, z: f64) -> f64; + fn get_densities(&self) -> &Vec; + fn get_concentrations(&self) -> &Vec; + fn get_electronic_stopping_correction_factor(&self) -> f64; } impl Mesh2D { @@ -117,7 +121,7 @@ impl Mesh2D { } } } -impl Mesh for Mesh2D { +impl Geometry for Mesh2D { fn get_energy_barrier_thickness(&self) -> f64 { self.energy_barrier_thickness } @@ -217,7 +221,7 @@ impl Cell2D { } } -impl MeshElement for Cell2D { +impl GeometryElement for Cell2D { /// Determines whether this cell contains the point (x, y). fn contains(&self, x: f64, y: f64, z: f64) -> bool { self.triangle.contains(x, y) @@ -227,6 +231,18 @@ impl MeshElement for Cell2D { fn distance_to(&self, x: f64, y: f64, z: f64) -> f64 { self.triangle.distance_to(x, y) } + + fn get_densities(&self) -> &Vec { + &self.densities + } + + fn get_concentrations(&self) -> &Vec { + &self.concentrations + } + + fn get_electronic_stopping_correction_factor(&self) -> f64 { + self.electronic_stopping_correction_factor + } } /// A triangle in 2D, with points (x1, y1), (x2, y2), (x3, y3), and the three line segments bewtween them. From 69530d7db1b90a49a1276ae6e1d7c8863a24700c Mon Sep 17 00:00:00 2001 From: Jon Drobny Date: Sat, 27 Feb 2021 20:01:26 -0600 Subject: [PATCH 05/14] Major rework of all geometry components. This version compiles, reproduces correct distributions, reflection coefficients, and sputtering yield. --- src/bca.rs | 20 +++++------ src/material.rs | 91 +++++++++++++++++++------------------------------ src/mesh.rs | 29 +++++++++------- 3 files changed, 61 insertions(+), 79 deletions(-) diff --git a/src/bca.rs b/src/bca.rs index ee3cd20..05bf502 100644 --- a/src/bca.rs +++ b/src/bca.rs @@ -1,5 +1,4 @@ use super::*; -use geo::{Closest}; #[cfg(any(feature = "cpr_rootfinder_openblas", feature = "cpr_rootfinder_netlib", feature = "cpr_rootfinder_intel_mkl"))] use rcpr::chebyshev::*; @@ -154,18 +153,15 @@ pub fn single_ion_bca(particle: particle::Particle, material: &mate let reduced_energy: f64 = LINDHARD_REDUCED_ENERGY_PREFACTOR*a*Mb/(Ma+Mb)/Za/Zb*E; let estimated_range_of_recoils = (reduced_energy.powf(0.3) + 0.1).powf(3.)/n/a/a; - //TODO: CHange this 2D - if let Closest::SinglePoint(p2) = material.closest_point(particle_2.pos.x, particle_2.pos.y, particle_2.pos.z) { - let dx = p2.x() - particle_2.pos.x; - let dy = p2.y() - particle_2.pos.y; - let distance_to_surface = (dx*dx + dy*dy).sqrt(); - - if (distance_to_surface < estimated_range_of_recoils) & (particle_2.E > particle_2.Ec) { - particles.push(particle_2); - } - } else { - panic!("Numerical error: geometry algorithm failed to find distance from particle to surface.") + let (x2, y2, z2) = material.closest_point(particle_2.pos.x, particle_2.pos.y, particle_2.pos.z); + let dx = x2 - particle_2.pos.x; + let dy = y2 - particle_2.pos.y; + let distance_to_surface = (dx*dx + dy*dy).sqrt(); + + if (distance_to_surface < estimated_range_of_recoils) & (particle_2.E > particle_2.Ec) { + particles.push(particle_2); } + //If transferred energy > cutoff energy, add recoil to particle vector } else if options.track_recoils & (particle_2.E > particle_2.Ec) { particles.push(particle_2); diff --git a/src/material.rs b/src/material.rs index c9c0e85..ded5220 100644 --- a/src/material.rs +++ b/src/material.rs @@ -1,8 +1,4 @@ use super::*; -use geo::algorithm::contains::Contains; -use geo::algorithm::closest_point::ClosestPoint; -use geo::{point, Closest}; - /// Holds material input parameters from [material_params]. #[derive(Deserialize)] @@ -73,7 +69,7 @@ impl Material { let total_number_density: f64 = if self.inside(x, y, z) { self.geometry.get_total_density(x, y, z) } else { - self.geometry.nearest_cell_to(x, y, z).densities.iter().sum::() + self.geometry.nearest_to(x, y, z).densities.iter().sum::() }; return self.geometry.get_densities(x, y, z).iter().map(|&i| i / total_number_density).collect(); @@ -102,7 +98,7 @@ impl Material { if self.inside(x, y, z) { return self.geometry.get_ck(x, y, z); } else { - return self.geometry.nearest_cell_to(x, y, z).electronic_stopping_correction_factor; + return self.geometry.nearest_to(x, y, z).electronic_stopping_correction_factor; } } @@ -116,7 +112,7 @@ impl Material { if self.inside(x, y, z) { return self.geometry.get_densities(x, y, z).iter().sum::(); } else { - return self.geometry.nearest_cell_to(x, y, z).densities.iter().sum::(); + return self.geometry.nearest_to(x, y, z).densities.iter().sum::(); } //return self.n.iter().sum::(); } @@ -126,7 +122,7 @@ impl Material { if self.inside(x, y, z) { return &self.geometry.get_densities(x, y, z); } else { - return &self.geometry.nearest_cell_to(x, y, z).densities; + return &self.geometry.nearest_to(x, y, z).densities; } //return self.n.iter().sum::(); } @@ -137,7 +133,7 @@ impl Material { if self.geometry.inside(x, y, z) { true } else { - let nearest_cell = self.geometry.nearest_cell_to(x, y, z); + let nearest_cell = self.geometry.nearest_to(x, y, z); let distance = nearest_cell.distance_to(x, y, z); distance < self.geometry.get_energy_barrier_thickness() } @@ -147,23 +143,12 @@ impl Material { /// Determines whether a point (x, y) is inside the simulation boundary. pub fn inside_simulation_boundary(&self, x:f64, y: f64, z: f64) -> bool { - let p = point!(x: x, y: y); - return self.geometry.get_simulation_boundary().contains(&p); + return self.geometry.inside_simulation_boundary(x, y, z); } - //TODO: CHANGE THIS TO REMOVE IMPLIED 2D /// Finds the closest point on the material boundary to the point (x, y). - pub fn closest_point(&self, x: f64, y: f64, z: f64) -> Closest { - let p = point!(x: x, y: y); - return self.geometry.get_boundary().closest_point(&p); - //return self.surface.closest_point(&p) - } - - //TODO: CHANGE THIS TO REMOVE IMPLIED 2D - /// Finds the closest point on the simulation boundary to the point (x, y). - pub fn closest_point_on_simulation_surface(&self, x: f64, y: f64, z: f64) -> Closest { - let p = point!(x: x, y: y); - return self.geometry.get_simulation_boundary().closest_point(&p) + pub fn closest_point(&self, x: f64, y: f64, z: f64) -> (f64, f64, f64) { + self.geometry.closest_point(x, y, z) } /// Finds the average, concentration-weighted atomic number, Z_effective, of the triangle that contains or is nearest to (x, y). @@ -331,10 +316,13 @@ pub fn surface_binding_energy(particle_1: &mut particle::Particle, if entering { if particle_1.backreflected { particle_1.backreflected = false; - } else if let Closest::SinglePoint(p2) = material.closest_point(x, y, z) { - let dx = p2.x() - x; - let dy = p2.y() - y; - let mag = (dx*dx + dy*dy).sqrt(); + } else { + let (x2, y2, z2) = material.closest_point(x, y, z); + + let dx = x2 - x; + let dy = y2 - y; + let dz = z2 - z; + let mag = (dx*dx + dy*dy + dz*dz).sqrt(); let costheta = dx*cosx/mag + dy*cosy/mag; particle_1.E += Es; @@ -342,47 +330,40 @@ pub fn surface_binding_energy(particle_1: &mut particle::Particle, //Surface refraction via Snell's law let delta_theta = particle::refraction_angle(costheta, E, E + Es); particle::rotate_particle(particle_1, delta_theta, 0.); - - //particle_1.add_trajectory(); - } else { - panic!("Numerical error: surface boundary algorithm encountered an error. Check geometry."); } } if leaving { - if let Closest::SinglePoint(p2) = material.closest_point(x, y, z) { - - let dx = p2.x() - x; - let dy = p2.y() - y; - let mag = (dx*dx + dy*dy).sqrt(); - let costheta = dx*cosx/mag + dy*cosy/mag; - let leaving_energy = E*costheta*costheta; + let (x2, y2, z2) = material.closest_point(x, y, z); + let dx = x2 - x; + let dy = y2 - y; + let dz = z2 - z; + let mag = (dx*dx + dy*dy + dz*dz).sqrt(); + let costheta = dx*cosx/mag + dy*cosy/mag; + let leaving_energy = E*costheta*costheta; - if costheta < 0. { - if leaving_energy > Es { + if costheta < 0. { + if leaving_energy > Es { - //particle_1.left = true; - particle_1.E += -Es; + //particle_1.left = true; + particle_1.E += -Es; - //Surface refraction via Snell's law - let delta_theta = particle::refraction_angle(costheta, E, E - Es); - particle::rotate_particle(particle_1, delta_theta, 0.); + //Surface refraction via Snell's law + let delta_theta = particle::refraction_angle(costheta, E, E - Es); + particle::rotate_particle(particle_1, delta_theta, 0.); - //particle_1.add_trajectory(); + //particle_1.add_trajectory(); - } else { + } else { - //Specular reflection at local surface normal - particle_1.dir.x = -2.*(costheta)*dx/mag + cosx; - particle_1.dir.y = -2.*(costheta)*dy/mag + cosy; + //Specular reflection at local surface normal + particle_1.dir.x = -2.*(costheta)*dx/mag + cosx; + particle_1.dir.y = -2.*(costheta)*dy/mag + cosy; - particle_1.backreflected = true; - //particle_1.add_trajectory(); - } + particle_1.backreflected = true; + //particle_1.add_trajectory(); } - } else { - panic!("Numerical error: surface boundary algorithm encountered an error. Check geometry."); } } } diff --git a/src/mesh.rs b/src/mesh.rs index 6bb338e..fb21a9e 100644 --- a/src/mesh.rs +++ b/src/mesh.rs @@ -1,7 +1,7 @@ use super::*; use geo::algorithm::contains::Contains; -use geo::{Polygon, LineString, Point}; - +use geo::{Polygon, LineString, Point, point, Closest}; +use geo::algorithm::closest_point::ClosestPoint; /// Object that contains raw mesh input data. #[derive(Deserialize)] @@ -29,11 +29,11 @@ pub trait Geometry { fn get_ck(&self, x: f64, y: f64, z: f64) -> f64; fn get_total_density(&self, x: f64, y: f64, z: f64) -> f64; fn get_concentrations(&self, x: f64, y: f64, z: f64) -> &Vec; - fn nearest_cell_to(&self, x: f64, y: f64, z: f64) -> &Cell2D; + fn nearest_to(&self, x: f64, y: f64, z: f64) -> &Cell2D; fn inside(&self, x: f64, y: f64, z: f64) -> bool; + fn inside_simulation_boundary(&self, x: f64, y: f64, z: f64) -> bool; fn get_energy_barrier_thickness(&self) -> f64; - fn get_simulation_boundary(&self) -> &Polygon; - fn get_boundary(&self) -> &Polygon; + fn closest_point(&self, x: f64, y: f64, z: f64) -> (f64, f64, f64); } pub trait GeometryElement { @@ -121,17 +121,23 @@ impl Mesh2D { } } } + impl Geometry for Mesh2D { fn get_energy_barrier_thickness(&self) -> f64 { self.energy_barrier_thickness } - fn get_simulation_boundary(&self) -> &Polygon { - &self.simulation_boundary + fn inside_simulation_boundary(&self, x: f64, y: f64, z: f64) -> bool { + self.simulation_boundary.contains(&point!(x: x, y: y)) } - fn get_boundary(&self) -> &Polygon { - &self.boundary + fn closest_point(&self, x: f64, y: f64, z: f64) -> (f64, f64, f64) { + if let Closest::SinglePoint(p) = self.boundary.closest_point(&point!(x: x, y: y)) { + (p.x(), p.y(), z) + } else { + panic!("Geometry error: closest point routine failed to find single closest point to ({}, {}, {}).", x, y, z); + } + } /// Find the number densities of the triangle that contains or is nearest to (x, y). @@ -174,7 +180,7 @@ impl Geometry for Mesh2D { } panic!("Geometry error: method inside() is returning true for points outside all cells. Check boundary points.") } else { - return &self.nearest_cell_to(x, y, z).concentrations; + return &self.nearest_to(x, y, z).concentrations; } } @@ -184,7 +190,7 @@ impl Geometry for Mesh2D { } /// Finds the cell that is nearest to (x, y). - fn nearest_cell_to(&self, x: f64, y: f64, z: f64) -> &Cell2D { + fn nearest_to(&self, x: f64, y: f64, z: f64) -> &Cell2D { let mut min_distance: f64 = std::f64::MAX; let mut index: usize = 0; @@ -201,7 +207,6 @@ impl Geometry for Mesh2D { } } - /// A mesh cell that contains a triangle and the local number densities and concentrations. pub struct Cell2D { triangle: Triangle2D, From 4c390d7bb6136d36bac098d23d7e01f1cb70e0d2 Mon Sep 17 00:00:00 2001 From: Jon Drobny Date: Sat, 27 Feb 2021 20:30:48 -0600 Subject: [PATCH 06/14] Removed some explicit 2D dependences in material.rs --- src/bca.rs | 5 +++-- src/material.rs | 22 ++++++------------- src/mesh.rs | 56 +++++++++++++++++++++++++++++++++---------------- 3 files changed, 47 insertions(+), 36 deletions(-) diff --git a/src/bca.rs b/src/bca.rs index 05bf502..9c48183 100644 --- a/src/bca.rs +++ b/src/bca.rs @@ -156,7 +156,8 @@ pub fn single_ion_bca(particle: particle::Particle, material: &mate let (x2, y2, z2) = material.closest_point(particle_2.pos.x, particle_2.pos.y, particle_2.pos.z); let dx = x2 - particle_2.pos.x; let dy = y2 - particle_2.pos.y; - let distance_to_surface = (dx*dx + dy*dy).sqrt(); + let dz = z2 - particle_2.pos.z; + let distance_to_surface = (dx*dx + dy*dy + dz*dz).sqrt(); if (distance_to_surface < estimated_range_of_recoils) & (particle_2.E > particle_2.Ec) { particles.push(particle_2); @@ -179,7 +180,7 @@ pub fn single_ion_bca(particle: particle::Particle, material: &mate strong_collision_index, &options); //Check boundary conditions on leaving and stopping - material::boundary_condition_2D_planar(&mut particle_1, &material); + material::boundary_condition_planar(&mut particle_1, &material); //Set particle index to topmost particle particle_index = particles.len(); diff --git a/src/material.rs b/src/material.rs index ded5220..f5a2a2c 100644 --- a/src/material.rs +++ b/src/material.rs @@ -69,7 +69,7 @@ impl Material { let total_number_density: f64 = if self.inside(x, y, z) { self.geometry.get_total_density(x, y, z) } else { - self.geometry.nearest_to(x, y, z).densities.iter().sum::() + self.geometry.get_densities_nearest_to(x, y, z).iter().sum::() }; return self.geometry.get_densities(x, y, z).iter().map(|&i| i / total_number_density).collect(); @@ -98,7 +98,7 @@ impl Material { if self.inside(x, y, z) { return self.geometry.get_ck(x, y, z); } else { - return self.geometry.nearest_to(x, y, z).electronic_stopping_correction_factor; + return self.geometry.get_ck_nearest_to(x, y, z); } } @@ -112,7 +112,7 @@ impl Material { if self.inside(x, y, z) { return self.geometry.get_densities(x, y, z).iter().sum::(); } else { - return self.geometry.nearest_to(x, y, z).densities.iter().sum::(); + return self.geometry.get_densities_nearest_to(x, y, z).iter().sum::(); } //return self.n.iter().sum::(); } @@ -122,23 +122,14 @@ impl Material { if self.inside(x, y, z) { return &self.geometry.get_densities(x, y, z); } else { - return &self.geometry.nearest_to(x, y, z).densities; + return &self.geometry.get_densities_nearest_to(x, y, z); } //return self.n.iter().sum::(); } /// Determines whether a point (x, y) is inside the energy barrier of the material. pub fn inside_energy_barrier(&self, x: f64, y: f64, z: f64) -> bool { - - if self.geometry.inside(x, y, z) { - true - } else { - let nearest_cell = self.geometry.nearest_to(x, y, z); - let distance = nearest_cell.distance_to(x, y, z); - distance < self.geometry.get_energy_barrier_thickness() - } - //let p = point!(x: x, y: y); - //return self.energy_surface.contains(&p); + self.geometry.inside_energy_barrier(x, y, z) } /// Determines whether a point (x, y) is inside the simulation boundary. @@ -312,7 +303,6 @@ pub fn surface_binding_energy(particle_1: &mut particle::Particle, let leaving = !inside_now & inside_old; let entering = inside_now & !inside_old; - //TODO: Remove explicit 2D dependence if entering { if particle_1.backreflected { particle_1.backreflected = false; @@ -369,7 +359,7 @@ pub fn surface_binding_energy(particle_1: &mut particle::Particle, } /// Apply the boundary conditions of a material on a particle, including stopping, leaving, and reflection/refraction by/through the surface binding potential. -pub fn boundary_condition_2D_planar(particle_1: &mut particle::Particle, material: &material::Material) { +pub fn boundary_condition_planar(particle_1: &mut particle::Particle, material: &material::Material) { let x = particle_1.pos.x; let y = particle_1.pos.y; let z = particle_1.pos.z; diff --git a/src/mesh.rs b/src/mesh.rs index fb21a9e..1e41da3 100644 --- a/src/mesh.rs +++ b/src/mesh.rs @@ -29,9 +29,11 @@ pub trait Geometry { fn get_ck(&self, x: f64, y: f64, z: f64) -> f64; fn get_total_density(&self, x: f64, y: f64, z: f64) -> f64; fn get_concentrations(&self, x: f64, y: f64, z: f64) -> &Vec; - fn nearest_to(&self, x: f64, y: f64, z: f64) -> &Cell2D; + fn get_densities_nearest_to(&self, x: f64, y: f64, z: f64) -> &Vec; + fn get_ck_nearest_to(&self, x: f64, y: f64, z: f64) -> f64; fn inside(&self, x: f64, y: f64, z: f64) -> bool; fn inside_simulation_boundary(&self, x: f64, y: f64, z: f64) -> bool; + fn inside_energy_barrier(&self, x: f64, y: f64, z: f64) -> bool; fn get_energy_barrier_thickness(&self) -> f64; fn closest_point(&self, x: f64, y: f64, z: f64) -> (f64, f64, f64); } @@ -120,9 +122,44 @@ impl Mesh2D { energy_barrier_thickness: energy_barrier_thickness, } } + + /// Finds the cell that is nearest to (x, y). + fn nearest_to(&self, x: f64, y: f64, z: f64) -> &Cell2D { + + let mut min_distance: f64 = std::f64::MAX; + let mut index: usize = 0; + + for (cell_index, cell) in self.mesh.iter().enumerate() { + let distance_to = cell.distance_to(x, y, z); + if distance_to < min_distance { + min_distance = distance_to; + index = cell_index; + } + } + + return &self.mesh[index]; + } } impl Geometry for Mesh2D { + fn inside_energy_barrier(&self, x: f64, y: f64, z: f64) -> bool { + if self.inside(x, y, z) { + true + } else { + let nearest_cell = self.nearest_to(x, y, z); + let distance = nearest_cell.distance_to(x, y, z); + distance < self.get_energy_barrier_thickness() + } + } + + fn get_ck_nearest_to(&self, x: f64, y: f64, z: f64) -> f64 { + self.nearest_to(x, y, z).get_electronic_stopping_correction_factor() + } + + fn get_densities_nearest_to(&self, x: f64, y: f64, z: f64) -> &Vec { + self.nearest_to(x, y, z).get_densities() + } + fn get_energy_barrier_thickness(&self) -> f64 { self.energy_barrier_thickness } @@ -188,23 +225,6 @@ impl Geometry for Mesh2D { fn inside(&self, x: f64, y: f64, z: f64) -> bool { self.boundary.contains(&Point::new(x, y)) } - - /// Finds the cell that is nearest to (x, y). - fn nearest_to(&self, x: f64, y: f64, z: f64) -> &Cell2D { - - let mut min_distance: f64 = std::f64::MAX; - let mut index: usize = 0; - - for (cell_index, cell) in self.mesh.iter().enumerate() { - let distance_to = cell.distance_to(x, y, z); - if distance_to < min_distance { - min_distance = distance_to; - index = cell_index; - } - } - - return &self.mesh[index]; - } } /// A mesh cell that contains a triangle and the local number densities and concentrations. From 930ffdf9f2504dcc06f844f419fdb384caa55acc Mon Sep 17 00:00:00 2001 From: Jon Drobny Date: Sat, 27 Feb 2021 22:10:41 -0600 Subject: [PATCH 07/14] First attempt to extend specular reflection and surface refraction to 3D. Removed old commented out code throughout, and added a couple why-comments. --- src/bca.rs | 3 +-- src/material.rs | 20 +++++--------------- src/particle.rs | 1 - 3 files changed, 6 insertions(+), 18 deletions(-) diff --git a/src/bca.rs b/src/bca.rs index 9c48183..96654c9 100644 --- a/src/bca.rs +++ b/src/bca.rs @@ -204,7 +204,6 @@ pub fn determine_mfp_phi_impact_parameter(particle_1: &mut particle let mut binary_collision_geometries = Vec::with_capacity(options.weak_collision_order + 1); //Each weak collision gets its own aziumuthal angle in annuli around collision point - //azimuthal angle randomly selected (0..2pi) for k in 0..options.weak_collision_order + 1 { phis_azimuthal.push(2.*PI*rand::random::()); } @@ -242,7 +241,7 @@ pub fn determine_mfp_phi_impact_parameter(particle_1: &mut particle pmax = (1./(material.total_number_density(x, y, z)*PI*ffp)).sqrt() } - //If free-flight-path less than the interatomic spacing, revert to solid model + //If free-flight-path less than the interatomic spacing, revert to amorphous solid model //Mentioned in Eckstein 1991, Ziegler, Biersack, and Ziegler 2008 (SRIM textbook 7-8) if ffp < mfp { diff --git a/src/material.rs b/src/material.rs index f5a2a2c..c9d3286 100644 --- a/src/material.rs +++ b/src/material.rs @@ -114,7 +114,6 @@ impl Material { } else { return self.geometry.get_densities_nearest_to(x, y, z).iter().sum::(); } - //return self.n.iter().sum::(); } /// Lists number density of each species of triangle that contains or is nearest to (x, y). @@ -124,7 +123,6 @@ impl Material { } else { return &self.geometry.get_densities_nearest_to(x, y, z); } - //return self.n.iter().sum::(); } /// Determines whether a point (x, y) is inside the energy barrier of the material. @@ -146,14 +144,12 @@ impl Material { pub fn average_Z(&self, x: f64, y: f64, z: f64) -> f64 { let concentrations = self.geometry.get_concentrations(x, y, z); return self.Z.iter().zip(concentrations).map(|(charge, concentration)| charge*concentration).collect::>().iter().sum(); - //return self.Z.iter().sum::()/self.Z.len() as f64; } /// Finds the average, concentration-weighted atomic mass, m_effective, of the triangle that contains or is nearest to (x, y). pub fn average_mass(&self, x: f64, y: f64, z: f64) -> f64 { let concentrations = self.geometry.get_concentrations(x, y, z); return self.m.iter().zip(concentrations).map(|(mass, concentration)| mass*concentration).collect::>().iter().sum(); - //return self.m.iter().sum::()/self.m.len() as f64; } /// Finds the average, concentration-weighted bulk binding energy of the triangle that contains or is nearest to (x, y). @@ -161,7 +157,6 @@ impl Material { //returns average bulk binding energy let concentrations = self.geometry.get_concentrations(x, y, z); return self.Eb.iter().zip(concentrations).map(|(bulk_binding_energy, concentration)| bulk_binding_energy*concentration).collect::>().iter().sum(); - //return self.Eb.iter().sum::()/self.Eb.len() as f64; } /// Finds the concentration-dependent surface binding energy of the triangle that contains or is nearest to (x, y). @@ -236,12 +231,11 @@ impl Material { let ck = self.electronic_stopping_correction_factor(x, y, z); for (n, Zb) in self.number_densities(x, y, z).iter().zip(&self.Z) { - //let n = self.number_density(particle_1.pos.x, particle_1.pos.y); - //let Zb = self.Z_eff(particle_1.pos.x, particle_1.pos.y); let beta = (1. - (1. + E/Ma/C.powf(2.)).powf(-2.)).sqrt(); let v = beta*C; + // This term is an empirical fit to the mean ionization potential let I0 = match *Zb < 13. { true => 12. + 7./Zb, false => 9.76 + 58.5*Zb.powf(-1.19), @@ -308,16 +302,15 @@ pub fn surface_binding_energy(particle_1: &mut particle::Particle, particle_1.backreflected = false; } else { let (x2, y2, z2) = material.closest_point(x, y, z); - let dx = x2 - x; let dy = y2 - y; let dz = z2 - z; let mag = (dx*dx + dy*dy + dz*dz).sqrt(); - let costheta = dx*cosx/mag + dy*cosy/mag; + let costheta = dx*cosx/mag + dy*cosy/mag + dz*cosz/mag; particle_1.E += Es; - //Surface refraction via Snell's law + //Surface refraction through energy barrier via Snell's law let delta_theta = particle::refraction_angle(costheta, E, E + Es); particle::rotate_particle(particle_1, delta_theta, 0.); } @@ -330,29 +323,26 @@ pub fn surface_binding_energy(particle_1: &mut particle::Particle, let dy = y2 - y; let dz = z2 - z; let mag = (dx*dx + dy*dy + dz*dz).sqrt(); - let costheta = dx*cosx/mag + dy*cosy/mag; + let costheta = dx*cosx/mag + dy*cosy/mag + dz*cosz/mag; let leaving_energy = E*costheta*costheta; if costheta < 0. { if leaving_energy > Es { - //particle_1.left = true; particle_1.E += -Es; //Surface refraction via Snell's law let delta_theta = particle::refraction_angle(costheta, E, E - Es); particle::rotate_particle(particle_1, delta_theta, 0.); - //particle_1.add_trajectory(); - } else { //Specular reflection at local surface normal particle_1.dir.x = -2.*(costheta)*dx/mag + cosx; particle_1.dir.y = -2.*(costheta)*dy/mag + cosy; + particle_1.dir.z = -2.*(costheta)*dz/mag + cosz; particle_1.backreflected = true; - //particle_1.add_trajectory(); } } } diff --git a/src/particle.rs b/src/particle.rs index 95ef3f4..62b3e39 100644 --- a/src/particle.rs +++ b/src/particle.rs @@ -210,7 +210,6 @@ pub fn particle_advance(particle_1: &mut particle::Particle, mfp: f64, asymptoti /// Calcualte the refraction angle based on the surface binding energy of the material. pub fn refraction_angle(costheta: f64, energy_old: f64, energy_new: f64) -> f64 { - //println!("energy_old: {} energy_new: {} costheta: {}", energy_old/EV, energy_new/EV, costheta); let costheta = if costheta.abs() > 1. {costheta.signum()} else {costheta}; let sintheta0 = (1. - costheta*costheta).sqrt(); let sintheta1 = sintheta0*(energy_old/energy_new).sqrt(); From 09777f3c18b68e45064ade7725143b8671314581 Mon Sep 17 00:00:00 2001 From: Jon Drobny Date: Sun, 28 Feb 2021 11:19:09 -0600 Subject: [PATCH 08/14] Additional under-the-hood changes to prepare for geometry abstraction. This was very difficult, and it's not in an amazing state, but it works and it compiles. The worst offender is not knowing how to share associated types between traits, so to make the material I have to clone and use the entire input file to avoid a circular reference between the InputFile and Geometry traits. --- examples/boron_nitride.toml | 2 +- examples/layered_geometry.toml | 4 +- examples/multiple_interaction_potentials.toml | 2 +- src/input.rs | 79 ++++++- src/main.rs | 4 +- src/material.rs | 11 +- src/mesh.rs | 201 +++++++++++++----- src/particle.rs | 2 +- src/tests.rs | 96 ++++++++- 9 files changed, 325 insertions(+), 76 deletions(-) diff --git a/examples/boron_nitride.toml b/examples/boron_nitride.toml index 8359226..047951e 100644 --- a/examples/boron_nitride.toml +++ b/examples/boron_nitride.toml @@ -44,7 +44,7 @@ dir = [ [ 0.9999999999984769, 1.7453292519934434e-6, 0.0,], [ 0.9999999999984769 interaction_index = [0, 0, 0, 0, 0, 0, 0, 0, 0, 0] particle_input_filename="" -[mesh_2d_input] +[geometry_input] length_unit = "MICRON" energy_barrier_thickness = 1.7645653881793786e-4 triangles = [ [ 0.0, 0.0, 0.0012558103905862675, 0.0, 0.02, 0.019960534568565433,], [ 0.0, 0.0012558103905862675, 0.0025066646712860853, 0.0, 0.019960534568565433, 0.019842294026289557,], [ 0.0, 0.0025066646712860853, 0.003747626291714493, 0.0, 0.019842294026289557, 0.019645745014573775,], [ 0.0, 0.003747626291714493, 0.004973797743297096, 0.0, 0.019645745014573775, 0.01937166322257262,], [ 0.0, 0.004973797743297096, 0.0061803398874989484, 0.0, 0.01937166322257262, 0.019021130325903073,], [ 0.0, 0.0061803398874989484, 0.00736249105369356, 0.0, 0.019021130325903073, 0.018595529717765027,], [ 0.0, 0.00736249105369356, 0.008515585831301454, 0.0, 0.018595529717765027, 0.01809654104932039,], [ 0.0, 0.008515585831301454, 0.009635073482034308, 0.0, 0.01809654104932039, 0.01752613360087727,], [ 0.0, 0.009635073482034308, 0.010716535899579934, 0.0, 0.01752613360087727, 0.0168865585100403,], [ 0.0, 0.010716535899579934, 0.011755705045849463, 0.0, 0.0168865585100403, 0.016180339887498948,], [ 0.0, 0.011755705045849463, 0.012748479794973795, 0.0, 0.016180339887498948, 0.015410264855515783,], [ 0.0, 0.012748479794973795, 0.013690942118573775, 0.0, 0.015410264855515783, 0.014579372548428232,], [ 0.0, 0.013690942118573775, 0.014579372548428232, 0.0, 0.014579372548428232, 0.013690942118573773,], [ 0.0, 0.014579372548428232, 0.015410264855515785, 0.0, 0.013690942118573773, 0.012748479794973793,], [ 0.0, 0.015410264855515785, 0.016180339887498948, 0.0, 0.012748479794973793, 0.011755705045849461,], [ 0.0, 0.016180339887498948, 0.0168865585100403, 0.0, 0.011755705045849461, 0.01071653589957993,], [ 0.0, 0.0168865585100403, 0.017526133600877274, 0.0, 0.01071653589957993, 0.009635073482034304,], [ 0.0, 0.017526133600877274, 0.018096541049320392, 0.0, 0.009635073482034304, 0.008515585831301454,], [ 0.0, 0.018096541049320392, 0.01859552971776503, 0.0, 0.008515585831301454, 0.007362491053693557,], [ 0.0, 0.01859552971776503, 0.019021130325903073, 0.0, 0.007362491053693557, 0.006180339887498949,], [ 0.0, 0.019021130325903073, 0.01937166322257262, 0.0, 0.006180339887498949, 0.004973797743297095,], [ 0.0, 0.01937166322257262, 0.019645745014573775, 0.0, 0.004973797743297095, 0.0037476262917144902,], [ 0.0, 0.019645745014573775, 0.019842294026289557, 0.0, 0.0037476262917144902, 0.0025066646712860853,], [ 0.0, 0.019842294026289557, 0.019960534568565433, 0.0, 0.0025066646712860853, 0.001255810390586266,], [ 0.0, 0.019960534568565433, 0.02, 0.0, 0.001255810390586266, -3.2162452993532734e-18,], [ 0.0, 0.02, 0.019960534568565433, 0.0, -3.2162452993532734e-18, -0.001255810390586268,], [ 0.0, 0.019960534568565433, 0.019842294026289557, 0.0, -0.001255810390586268, -0.0025066646712860875,], [ 0.0, 0.019842294026289557, 0.019645745014573772, 0.0, -0.0025066646712860875, -0.0037476262917144967,], [ 0.0, 0.019645745014573772, 0.01937166322257262, 0.0, -0.0037476262917144967, -0.004973797743297097,], [ 0.0, 0.01937166322257262, 0.019021130325903073, 0.0, -0.004973797743297097, -0.006180339887498951,], [ 0.0, 0.019021130325903073, 0.01859552971776503, 0.0, -0.006180339887498951, -0.00736249105369356,], [ 0.0, 0.01859552971776503, 0.01809654104932039, 0.0, -0.00736249105369356, -0.008515585831301454,], [ 0.0, 0.01809654104932039, 0.01752613360087727, 0.0, -0.008515585831301454, -0.00963507348203431,], [ 0.0, 0.01752613360087727, 0.0168865585100403, 0.0, -0.00963507348203431, -0.010716535899579938,], [ 0.0, 0.0168865585100403, 0.016180339887498948, 0.0, -0.010716535899579938, -0.011755705045849461,], [ 0.0, 0.016180339887498948, 0.015410264855515785, 0.0, -0.011755705045849461, -0.012748479794973795,], [ 0.0, 0.015410264855515785, 0.014579372548428228, 0.0, -0.012748479794973795, -0.013690942118573775,], [ 0.0, 0.014579372548428228, 0.01369094211857377, 0.0, -0.013690942118573775, -0.014579372548428234,], [ 0.0, 0.01369094211857377, 0.01274847979497379, 0.0, -0.014579372548428234, -0.015410264855515788,], [ 0.0, 0.01274847979497379, 0.011755705045849465, 0.0, -0.015410264855515788, -0.016180339887498948,], [ 0.0, 0.011755705045849465, 0.010716535899579934, 0.0, -0.016180339887498948, -0.0168865585100403,], [ 0.0, 0.010716535899579934, 0.009635073482034304, 0.0, -0.0168865585100403, -0.01752613360087727,], [ 0.0, 0.009635073482034304, 0.00851558583130145, 0.0, -0.01752613360087727, -0.018096541049320392,], [ 0.0, 0.00851558583130145, 0.007362491053693555, 0.0, -0.018096541049320392, -0.01859552971776503,], [ 0.0, 0.007362491053693555, 0.006180339887498942, 0.0, -0.01859552971776503, -0.019021130325903073,], [ 0.0, 0.006180339887498942, 0.004973797743297097, 0.0, -0.019021130325903073, -0.01937166322257262,], [ 0.0, 0.004973797743297097, 0.0037476262917144915, 0.0, -0.01937166322257262, -0.019645745014573775,], [ 0.0, 0.0037476262917144915, 0.002506664671286082, 0.0, -0.019645745014573775, -0.019842294026289557,], [ 0.0, 0.002506664671286082, 0.0012558103905862628, 0.0, -0.019842294026289557, -0.019960534568565433,], [ 0.0, 0.0012558103905862628, -6.432490598706547e-18, 0.0, -0.019960534568565433, -0.02,], [ 0.0, -6.432490598706547e-18, -0.0012558103905862669, 0.0, -0.02, -0.019960534568565433,], [ 0.0, -0.0012558103905862669, -0.0025066646712860858, 0.0, -0.019960534568565433, -0.019842294026289557,], [ 0.0, -0.0025066646712860858, -0.0037476262917144954, 0.0, -0.019842294026289557, -0.019645745014573772,], [ 0.0, -0.0037476262917144954, -0.0049737977432971, 0.0, -0.019645745014573772, -0.01937166322257262,], [ 0.0, -0.0049737977432971, -0.0061803398874989545, 0.0, -0.01937166322257262, -0.019021130325903073,], [ 0.0, -0.0061803398874989545, -0.007362491053693567, 0.0, -0.019021130325903073, -0.018595529717765024,], [ 0.0, -0.007362491053693567, -0.008515585831301454, 0.0, -0.018595529717765024, -0.01809654104932039,], [ 0.0, -0.008515585831301454, -0.009635073482034308, 0.0, -0.01809654104932039, -0.01752613360087727,], [ 0.0, -0.009635073482034308, -0.010716535899579936, 0.0, -0.01752613360087727, -0.0168865585100403,], [ 0.0, -0.010716535899579936, -0.011755705045849468, 0.0, -0.0168865585100403, -0.016180339887498944,], [ 0.0, -0.011755705045849468, -0.0127484797949738, 0.0, -0.016180339887498944, -0.015410264855515781,], [ 0.0, -0.0127484797949738, -0.013690942118573775, 0.0, -0.015410264855515781, -0.014579372548428232,], [ 0.0, -0.013690942118573775, -0.014579372548428232, 0.0, -0.014579372548428232, -0.013690942118573773,], [ 0.0, -0.014579372548428232, -0.015410264855515788, 0.0, -0.013690942118573773, -0.01274847979497379,], [ 0.0, -0.015410264855515788, -0.016180339887498948, 0.0, -0.01274847979497379, -0.011755705045849465,], [ 0.0, -0.016180339887498948, -0.016886558510040308, 0.0, -0.011755705045849465, -0.010716535899579927,], [ 0.0, -0.016886558510040308, -0.01752613360087727, 0.0, -0.010716535899579927, -0.009635073482034306,], [ 0.0, -0.01752613360087727, -0.018096541049320396, 0.0, -0.009635073482034306, -0.008515585831301443,], [ 0.0, -0.018096541049320396, -0.01859552971776503, 0.0, -0.008515585831301443, -0.007362491053693556,], [ 0.0, -0.01859552971776503, -0.019021130325903073, 0.0, -0.007362491053693556, -0.006180339887498951,], [ 0.0, -0.019021130325903073, -0.019371663222572624, 0.0, -0.006180339887498951, -0.004973797743297089,], [ 0.0, -0.019371663222572624, -0.019645745014573775, 0.0, -0.004973797743297089, -0.003747626291714493,], [ 0.0, -0.019645745014573775, -0.019842294026289557, 0.0, -0.003747626291714493, -0.0025066646712860745,], [ 0.0, -0.019842294026289557, -0.019960534568565433, 0.0, -0.0025066646712860745, -0.001255810390586264,], [ 0.0, -0.019960534568565433, -0.02, 0.0, -0.001255810390586264, -3.673940397442059e-18,], [ 0.0, -0.02, -0.019960534568565433, 0.0, -3.673940397442059e-18, 0.0012558103905862745,], [ 0.0, -0.019960534568565433, -0.019842294026289557, 0.0, 0.0012558103905862745, 0.0025066646712860845,], [ 0.0, -0.019842294026289557, -0.019645745014573772, 0.0, 0.0025066646712860845, 0.003747626291714503,], [ 0.0, -0.019645745014573772, -0.01937166322257262, 0.0, 0.003747626291714503, 0.004973797743297099,], [ 0.0, -0.01937166322257262, -0.019021130325903073, 0.0, 0.004973797743297099, 0.006180339887498945,], [ 0.0, -0.019021130325903073, -0.018595529717765024, 0.0, 0.006180339887498945, 0.007362491053693565,], [ 0.0, -0.018595529717765024, -0.018096541049320392, 0.0, 0.007362491053693565, 0.008515585831301452,], [ 0.0, -0.018096541049320392, -0.017526133600877267, 0.0, 0.008515585831301452, 0.009635073482034314,], [ 0.0, -0.017526133600877267, -0.0168865585100403, 0.0, 0.009635073482034314, 0.010716535899579936,], [ 0.0, -0.0168865585100403, -0.01618033988749894, 0.0, 0.010716535899579936, 0.011755705045849473,], [ 0.0, -0.01618033988749894, -0.015410264855515781, 0.0, 0.011755705045849473, 0.0127484797949738,], [ 0.0, -0.015410264855515781, -0.014579372548428232, 0.0, 0.0127484797949738, 0.013690942118573773,], [ 0.0, -0.014579372548428232, -0.013690942118573766, 0.0, 0.013690942118573773, 0.014579372548428239,], [ 0.0, -0.013690942118573766, -0.012748479794973793, 0.0, 0.014579372548428239, 0.015410264855515788,], [ 0.0, -0.012748479794973793, -0.011755705045849453, 0.0, 0.015410264855515788, 0.016180339887498955,], [ 0.0, -0.011755705045849453, -0.010716535899579927, 0.0, 0.016180339887498955, 0.016886558510040305,], [ 0.0, -0.010716535899579927, -0.009635073482034308, 0.0, 0.016886558510040305, 0.01752613360087727,], [ 0.0, -0.009635073482034308, -0.008515585831301445, 0.0, 0.01752613360087727, 0.018096541049320396,], [ 0.0, -0.008515585831301445, -0.007362491053693557, 0.0, 0.018096541049320396, 0.01859552971776503,], [ 0.0, -0.007362491053693557, -0.006180339887498935, 0.0, 0.01859552971776503, 0.019021130325903076,], [ 0.0, -0.006180339887498935, -0.00497379774329709, 0.0, 0.019021130325903076, 0.019371663222572624,], [ 0.0, -0.00497379774329709, -0.0037476262917144937, 0.0, 0.019371663222572624, 0.019645745014573775,], [ 0.0, -0.0037476262917144937, -0.002506664671286076, 0.0, 0.019645745014573775, 0.019842294026289557,], [ 0.0, -0.002506664671286076, -0.0012558103905862654, 0.0, 0.019842294026289557, 0.019960534568565433,], [ 0.0, -0.0012558103905862654, -4.898587196589413e-18, 0.0, 0.019960534568565433, 0.02,],] diff --git a/examples/layered_geometry.toml b/examples/layered_geometry.toml index 5f22a46..01da35b 100644 --- a/examples/layered_geometry.toml +++ b/examples/layered_geometry.toml @@ -45,10 +45,10 @@ Ec = [ 1.0,] Es = [ 0.0,] interaction_index = [0] pos = [ [ -1.7453292519934434e-8, 0.0, 0.0,],] -dir = [ [ 0.999, 0.0001, 0.0,],] +dir = [ [ 0.999, 0.001, 0.0,],] particle_input_filename = "" -[mesh_2d_input] +[geometry_input] length_unit = "MICRON" triangles = [ [ 0.0, 0.01, 0.0, 0.5, 0.5, -0.5,], [ 0.0, 0.01, 0.01, -0.5, 0.5, -0.5,], [ 0.01, 0.01, 0.04, -0.5, 0.5, -0.5,], [ 0.01, 0.04, 0.04, 0.5, 0.5, -0.5,], [ 0.04, 0.5, 0.04, 0.5, 0.5, -0.5,], [ 0.04, 0.5, 0.5, -0.5, 0.5, -0.5,],] densities = [ [ 3.0e+10, 6.0e+10, 0.0, 0.0,], [ 3.0e+10, 6.0e+10, 0.0, 0.0,], [ 0.0, 0.0, 6.026e+10, 0.0,], [ 0.0, 0.0, 6.026e+10, 0.0,], [ 0.0, 0.0, 0.0, 4.996e+10,], [ 0.0, 0.0, 0.0, 4.996e+10,],] diff --git a/examples/multiple_interaction_potentials.toml b/examples/multiple_interaction_potentials.toml index f5a49ea..f10fa32 100644 --- a/examples/multiple_interaction_potentials.toml +++ b/examples/multiple_interaction_potentials.toml @@ -42,7 +42,7 @@ pos = [ [ -1.7453292519934434e-8, 0.0, 0.0,], [ -1.7453292519934434e-8, 0.0, 0.0 dir = [ [ 0.9999999999984769, 1.7453292519934434e-6, 0.0,], [ 0.9999999999984769, 1.7453292519934434e-6, 0.0,],] particle_input_filename = -[mesh_2d_input] +[geometry_input] energy_barrier_thickness = 2.2E-4 length_unit = "MICRON" triangles = [ [ 0.0, 0.01, 0.0, 0.5, 0.5, -0.5,], [ 0.0, 0.01, 0.01, -0.5, 0.5, -0.5,], [ 0.01, 0.01, 0.04, -0.5, 0.5, -0.5,], [ 0.01, 0.04, 0.04, 0.5, 0.5, -0.5,], [ 0.04, 0.5, 0.04, 0.5, 0.5, -0.5,], [ 0.04, 0.5, 0.5, -0.5, 0.5, -0.5,],] diff --git a/src/input.rs b/src/input.rs index 8b1db8b..b821750 100644 --- a/src/input.rs +++ b/src/input.rs @@ -1,18 +1,76 @@ use super::*; +pub trait InputFile { + type GeometryInput; + fn new(string: &str) -> Self; + fn get_options(&self) -> &Options; + fn get_material_parameters(&self) -> &material::MaterialParameters; + fn get_particle_parameters(&self) -> &particle::ParticleParameters; + fn get_geometry_input(&self) -> &Self::GeometryInput; +} /// Rustbca's internal representation of an input file. #[derive(Deserialize)] -pub struct Input { +pub struct Input2D { + pub options: Options, + pub material_parameters: material::MaterialParameters, + pub particle_parameters: particle::ParticleParameters, + pub geometry_input: mesh::Mesh2DInput, +} + +impl InputFile for Input2D { + type GeometryInput = mesh::Mesh2DInput; + + fn new(string: &str) -> Input2D { + toml::from_str(string).expect("Could not parse TOML file.") + } + + fn get_options(&self) -> &Options{ + &self.options + } + fn get_material_parameters(&self) -> &material::MaterialParameters{ + &self.material_parameters + } + fn get_particle_parameters(&self) -> &particle::ParticleParameters{ + &self.particle_parameters + } + fn get_geometry_input(&self) -> &Self::GeometryInput{ + &self.geometry_input + } +} + +#[derive(Deserialize)] +pub struct Input0D { pub options: Options, pub material_parameters: material::MaterialParameters, pub particle_parameters: particle::ParticleParameters, - pub mesh_2d_input: mesh::Mesh2DInput, + pub geometry_input: mesh::Mesh0DInput, +} + +impl InputFile for Input0D { + type GeometryInput = mesh::Mesh0DInput; + + fn new(string: &str) -> Input0D { + toml::from_str(string).expect("Could not parse TOML file.") + } + + fn get_options(&self) -> &Options{ + &self.options + } + fn get_material_parameters(&self) -> &material::MaterialParameters{ + &self.material_parameters + } + fn get_particle_parameters(&self) ->& particle::ParticleParameters{ + &self.particle_parameters + } + fn get_geometry_input(&self) -> &Self::GeometryInput{ + &self.geometry_input + } } /// Rustbca's internal representation of the simulation-level options. #[cfg(not(feature = "distributions"))] -#[derive(Deserialize)] +#[derive(Deserialize, Clone)] pub struct Options { pub name: String, pub track_trajectories: bool, @@ -35,7 +93,7 @@ pub struct Options { } #[cfg(feature = "distributions")] -#[derive(Deserialize)] +#[derive(Deserialize, Clone)] pub struct Options { pub name: String, pub track_trajectories: bool, @@ -72,7 +130,8 @@ pub struct Options { pub z_num: usize, } -pub fn input() -> (Vec, material::Material, Options, OutputUnits){ +pub fn input() -> (Vec, material::Material, Options, OutputUnits) +where ::InputFileFormat: Deserialize<'static> + 'static { let args: Vec = env::args().collect(); @@ -91,13 +150,13 @@ pub fn input() -> (Vec, material::Material::InputFileFormat = InputFile::new(&input_toml); //Unpack toml information into structs - let material: material::Material = material::Material::::new(input.material_parameters, input.mesh_2d_input); - let options = input.options; - let particle_parameters = input.particle_parameters; + let material: material::Material = material::Material::::new(&input); + let options = (*input.get_options()).clone(); + let particle_parameters = (*input.get_particle_parameters()).clone(); //Ensure nonsensical threads/chunks options crash on input assert!(options.num_threads > 0, "Input error: num_threads must be greater than zero."); diff --git a/src/main.rs b/src/main.rs index a1c88d4..428d002 100644 --- a/src/main.rs +++ b/src/main.rs @@ -62,13 +62,13 @@ pub mod structs; pub use crate::enums::*; pub use crate::consts::*; pub use crate::structs::*; -pub use crate::input::{Input, Options}; +pub use crate::input::{Input2D, Input0D, Options, InputFile}; pub use crate::output::{OutputUnits}; pub use crate::mesh::{Geometry, GeometryElement}; fn main() { //Open and process input_file - let (particle_input_array, material, options, output_units) = input::input(); + let (particle_input_array, material, options, output_units): (Vec, material::Material, Options, OutputUnits) = input::input(); println!("Processing {} ions...", particle_input_array.len()); diff --git a/src/material.rs b/src/material.rs index c9d3286..694fdfd 100644 --- a/src/material.rs +++ b/src/material.rs @@ -28,8 +28,9 @@ pub struct Material { } impl Material { - /// Constructs a new material object from a material parameters object and a mesh_2d_input object. - pub fn new(material_parameters: MaterialParameters, mesh_2d_input: mesh::Mesh2DInput) -> Material { + pub fn new(raw_input: &::InputFileFormat) -> Material { + + let material_parameters = raw_input.get_material_parameters(); let energy_unit: f64 = match material_parameters.energy_unit.as_str() { "EV" => EV, @@ -53,13 +54,13 @@ impl Material { Material { m: material_parameters.m.iter().map(|&i| i*mass_unit).collect(), - Z: material_parameters.Z, + Z: material_parameters.Z.clone(), Eb: material_parameters.Eb.iter().map(|&i| i*energy_unit).collect(), Es: material_parameters.Es.iter().map(|&i| i*energy_unit).collect(), Ec: material_parameters.Ec.iter().map(|&i| i*energy_unit).collect(), - interaction_index: material_parameters.interaction_index, - geometry: Box::new(mesh::Mesh2D::new(mesh_2d_input)), + interaction_index: material_parameters.interaction_index.clone(), surface_binding_model: material_parameters.surface_binding_model, + geometry: Box::new(T::new(raw_input)), } } diff --git a/src/mesh.rs b/src/mesh.rs index 1e41da3..e585c0a 100644 --- a/src/mesh.rs +++ b/src/mesh.rs @@ -3,28 +3,14 @@ use geo::algorithm::contains::Contains; use geo::{Polygon, LineString, Point, point, Closest}; use geo::algorithm::closest_point::ClosestPoint; -/// Object that contains raw mesh input data. -#[derive(Deserialize)] -pub struct Mesh2DInput { - pub length_unit: String, - pub triangles: Vec<(f64, f64, f64, f64, f64, f64)>, - pub material_boundary_points: Vec<(f64, f64)>, - pub simulation_boundary_points: Vec<(f64, f64)>, - pub densities: Vec>, - pub electronic_stopping_correction_factors: Vec, - pub energy_barrier_thickness: f64, -} - -/// Triangular mesh for rustbca. -pub struct Mesh2D { - mesh: Vec, - pub boundary: Polygon, - pub simulation_boundary: Polygon, - pub energy_barrier_thickness: f64 -} -///Trait for a Geometry object +///Trait for a Geometry object - all forms of geometry must implement these traits to be used pub trait Geometry { + type GeometryInput; + type InputFileFormat: InputFile; + + fn new(input: &Self::InputFileFormat) -> Self; + fn get_densities(&self, x: f64, y: f64, z: f64) -> &Vec; fn get_ck(&self, x: f64, y: f64, z: f64) -> f64; fn get_total_density(&self, x: f64, y: f64, z: f64) -> f64; @@ -46,33 +32,167 @@ pub trait GeometryElement { fn get_electronic_stopping_correction_factor(&self) -> f64; } +#[derive(Deserialize)] +pub struct Mesh0DInput { + pub length_unit: String, + pub densities: Vec, + pub electronic_stopping_correction_factor: f64, + pub energy_barrier_thickness: f64, +} + +pub struct Mesh0D { + pub densities: Vec, + pub concentrations: Vec, + pub electronic_stopping_correction_factor: f64, + pub energy_barrier_thickness: f64, +} + +impl Geometry for Mesh0D { + + type GeometryInput = Mesh0DInput; + type InputFileFormat = Input0D; + + fn new(input: &Self::InputFileFormat) -> Mesh0D { + + let input = input.get_geometry_input(); + + let length_unit: f64 = match input.length_unit.as_str() { + "MICRON" => MICRON, + "CM" => CM, + "ANGSTROM" => ANGSTROM, + "NM" => NM, + "M" => 1., + _ => input.length_unit.parse() + .expect(format!( + "Input errror: could nor parse length unit {}. Use a valid float or one of ANGSTROM, NM, MICRON, CM, MM, M", + &input.length_unit.as_str() + ).as_str()), + }; + + let electronic_stopping_correction_factor = input.electronic_stopping_correction_factor; + let energy_barrier_thickness = input.energy_barrier_thickness/length_unit; + + let densities: Vec = input.densities.iter().map(|element| element/(length_unit).powf(3.)).collect(); + + let total_density: f64 = densities.iter().sum(); + let concentrations: Vec = densities.iter().map(|&density| density/total_density).collect::>(); + + Mesh0D { + densities, + concentrations, + electronic_stopping_correction_factor, + energy_barrier_thickness + } + } + + fn get_densities(&self, x: f64, y: f64, z: f64) -> &Vec { + &self.densities + } + fn get_ck(&self, x: f64, y: f64, z: f64) -> f64 { + self.electronic_stopping_correction_factor + } + fn get_total_density(&self, x: f64, y: f64, z: f64) -> f64{ + self.densities.iter().sum() + } + fn get_concentrations(&self, x: f64, y: f64, z: f64) -> &Vec { + &self.concentrations + } + fn get_densities_nearest_to(&self, x: f64, y: f64, z: f64) -> &Vec { + &self.densities + } + fn get_ck_nearest_to(&self, x: f64, y: f64, z: f64) -> f64 { + self.electronic_stopping_correction_factor + } + fn inside(&self, x: f64, y: f64, z: f64) -> bool { + x > 0.0 + } + fn inside_simulation_boundary(&self, x: f64, y: f64, z: f64) -> bool { + x > -2.*self.energy_barrier_thickness + } + fn inside_energy_barrier(&self, x: f64, y: f64, z: f64) -> bool { + x > -self.energy_barrier_thickness + } + + fn get_energy_barrier_thickness(&self) -> f64 { + self.energy_barrier_thickness + } + + fn closest_point(&self, x: f64, y: f64, z: f64) -> (f64, f64, f64) { + (0., y, z) + } +} + +/// Object that contains raw mesh input data. +#[derive(Deserialize)] +pub struct Mesh2DInput { + pub length_unit: String, + pub triangles: Vec<(f64, f64, f64, f64, f64, f64)>, + pub material_boundary_points: Vec<(f64, f64)>, + pub simulation_boundary_points: Vec<(f64, f64)>, + pub densities: Vec>, + pub electronic_stopping_correction_factors: Vec, + pub energy_barrier_thickness: f64, +} + +/// Triangular mesh for rustbca. +pub struct Mesh2D { + mesh: Vec, + pub boundary: Polygon, + pub simulation_boundary: Polygon, + pub energy_barrier_thickness: f64 +} + impl Mesh2D { - /// Constructor for Mesh2D object from mesh_2d_input. - pub fn new(mesh_2d_input: Mesh2DInput) -> Mesh2D { + /// Finds the cell that is nearest to (x, y). + fn nearest_to(&self, x: f64, y: f64, z: f64) -> &Cell2D { - let triangles = mesh_2d_input.triangles; - let material_boundary_points = mesh_2d_input.material_boundary_points; - let simulation_boundary_points = mesh_2d_input.simulation_boundary_points; - let electronic_stopping_correction_factors = mesh_2d_input.electronic_stopping_correction_factors; + let mut min_distance: f64 = std::f64::MAX; + let mut index: usize = 0; + + for (cell_index, cell) in self.mesh.iter().enumerate() { + let distance_to = cell.distance_to(x, y, z); + if distance_to < min_distance { + min_distance = distance_to; + index = cell_index; + } + } + + return &self.mesh[index]; + } +} + +impl Geometry for Mesh2D { + + type InputFileFormat = Input2D; + type GeometryInput = Mesh2DInput; + + /// Constructor for Mesh2D object from geometry_input. + fn new(input: &Self::InputFileFormat) -> Mesh2D { + + let geometry_input = input.get_geometry_input(); + let triangles = geometry_input.triangles.clone(); + let material_boundary_points = geometry_input.material_boundary_points.clone(); + let simulation_boundary_points = geometry_input.simulation_boundary_points.clone(); + let electronic_stopping_correction_factors = geometry_input.electronic_stopping_correction_factors.clone(); let n = triangles.len(); let mut cells: Vec = Vec::with_capacity(n); //Multiply all coordinates by value of geometry unit. - let length_unit: f64 = match mesh_2d_input.length_unit.as_str() { + let length_unit: f64 = match geometry_input.length_unit.as_str() { "MICRON" => MICRON, "CM" => CM, "ANGSTROM" => ANGSTROM, "NM" => NM, "M" => 1., - _ => mesh_2d_input.length_unit.parse() + _ => geometry_input.length_unit.parse() .expect(format!( - "Input errror: could nor parse length unit {}. Use a valid float or one of ANGSTROM, NM, MICRON, CM, MM, M", &mesh_2d_input.length_unit.as_str() + "Input errror: could nor parse length unit {}. Use a valid float or one of ANGSTROM, NM, MICRON, CM, MM, M", &geometry_input.length_unit.as_str() ).as_str()), }; - let densities: Vec> = mesh_2d_input.densities + let densities: Vec> = geometry_input.densities .iter() .map( |row| row.iter().map(|element| element/(length_unit).powf(3.)).collect() ).collect(); @@ -113,7 +233,7 @@ impl Mesh2D { } let simulation_boundary: Polygon = Polygon::new(LineString::from(simulation_boundary_points_converted), vec![]); - let energy_barrier_thickness = mesh_2d_input.energy_barrier_thickness*length_unit; + let energy_barrier_thickness = geometry_input.energy_barrier_thickness*length_unit; Mesh2D { mesh: cells, @@ -123,25 +243,6 @@ impl Mesh2D { } } - /// Finds the cell that is nearest to (x, y). - fn nearest_to(&self, x: f64, y: f64, z: f64) -> &Cell2D { - - let mut min_distance: f64 = std::f64::MAX; - let mut index: usize = 0; - - for (cell_index, cell) in self.mesh.iter().enumerate() { - let distance_to = cell.distance_to(x, y, z); - if distance_to < min_distance { - min_distance = distance_to; - index = cell_index; - } - } - - return &self.mesh[index]; - } -} - -impl Geometry for Mesh2D { fn inside_energy_barrier(&self, x: f64, y: f64, z: f64) -> bool { if self.inside(x, y, z) { true diff --git a/src/particle.rs b/src/particle.rs index 62b3e39..8387417 100644 --- a/src/particle.rs +++ b/src/particle.rs @@ -1,7 +1,7 @@ use super::*; /// Rustbca's internal representation of the particle_parameters input. -#[derive(Deserialize)] +#[derive(Deserialize, Clone)] pub struct ParticleParameters { pub particle_input_filename: String, pub length_unit: String, diff --git a/src/tests.rs b/src/tests.rs index 5db578b..ddf1fa5 100644 --- a/src/tests.rs +++ b/src/tests.rs @@ -32,7 +32,7 @@ fn test_surface_binding_energy_barrier() { let thickness: f64 = 1000.; let depth: f64 = 1000.; - let mesh_2d_input = mesh::Mesh2DInput { + let geometry_input = mesh::Mesh2DInput { length_unit: "ANGSTROM".to_string(), triangles: vec![(0., depth, 0., thickness/2., thickness/2., -thickness/2.), (depth, depth, 0., thickness/2., -thickness/2., -thickness/2.)], densities: vec![vec![0.03, 0.03], vec![0.03, 0.03]], @@ -42,7 +42,51 @@ fn test_surface_binding_energy_barrier() { electronic_stopping_correction_factors: vec![0.0, 0.0], }; - let material_1: material::Material = material::Material::::new(material_parameters, mesh_2d_input); + let options = Options { + name: "test".to_string(), + track_trajectories: false, + track_recoils: true, + track_recoil_trajectories: false, + write_buffer_size: 8000, + weak_collision_order: 0, + suppress_deep_recoils: false, + high_energy_free_flight_paths: true, + electronic_stopping_mode: ElectronicStoppingMode::INTERPOLATED, + mean_free_path_model: MeanFreePathModel::LIQUID, + interaction_potential: vec![vec![InteractionPotential::MOLIERE]], + scattering_integral: vec![vec![ScatteringIntegral::MENDENHALL_WELLER]], + num_threads: 1, + num_chunks: 1, + use_hdf5: false, + root_finder: vec![vec![Rootfinder::NEWTON{max_iterations: 100, tolerance: 1E-3}]], + track_displacements: false, + track_energy_losses: false, + }; + + let particle_parameters = particle::ParticleParameters { + particle_input_filename: "test".to_string(), + length_unit: "ANGSTROM".to_string(), + energy_unit: "EV".to_string(), + mass_unit: "AMU".to_string(), + N: vec![], + m: vec![], + Z: vec![], + E: vec![], + Ec: vec![], + Es: vec![], + pos: vec![], + dir: vec![], + interaction_index: vec![], + }; + + let raw_input = Input2D { + options, + geometry_input, + material_parameters, + particle_parameters, + }; + + let material_1: material::Material = material::Material::::new(&raw_input); particle_1.pos.x = 500.*ANGSTROM; particle_1.pos.y = 0.; @@ -215,7 +259,7 @@ fn test_momentum_conservation() { let thickness: f64 = 1000.; let depth: f64 = 1000.; - let mesh_2d_input = mesh::Mesh2DInput { + let geometry_input = mesh::Mesh2DInput { length_unit: "ANGSTROM".to_string(), triangles: vec![(0., depth, 0., thickness/2., thickness/2., -thickness/2.), (depth, depth, 0., thickness/2., -thickness/2., -thickness/2.)], densities: vec![vec![0.06306], vec![0.06306]], @@ -225,7 +269,51 @@ fn test_momentum_conservation() { energy_barrier_thickness: 0., }; - let material_1: material::Material = material::Material::::new(material_parameters, mesh_2d_input); + let options = Options { + name: "test".to_string(), + track_trajectories: false, + track_recoils: true, + track_recoil_trajectories: false, + write_buffer_size: 8000, + weak_collision_order: 0, + suppress_deep_recoils: false, + high_energy_free_flight_paths: true, + electronic_stopping_mode: ElectronicStoppingMode::INTERPOLATED, + mean_free_path_model: MeanFreePathModel::LIQUID, + interaction_potential: vec![vec![InteractionPotential::MOLIERE]], + scattering_integral: vec![vec![ScatteringIntegral::MENDENHALL_WELLER]], + num_threads: 1, + num_chunks: 1, + use_hdf5: false, + root_finder: vec![vec![Rootfinder::NEWTON{max_iterations: 100, tolerance: 1E-3}]], + track_displacements: false, + track_energy_losses: false, + }; + + let particle_parameters = particle::ParticleParameters { + particle_input_filename: "test".to_string(), + length_unit: "ANGSTROM".to_string(), + energy_unit: "EV".to_string(), + mass_unit: "AMU".to_string(), + N: vec![], + m: vec![], + Z: vec![], + E: vec![], + Ec: vec![], + Es: vec![], + pos: vec![], + dir: vec![], + interaction_index: vec![], + }; + + let raw_input = Input2D { + options, + geometry_input, + material_parameters, + particle_parameters, + }; + + let material_1: material::Material = material::Material::::new(&raw_input); for high_energy_free_flight_paths in vec![true, false] { for potential in vec![InteractionPotential::KR_C, InteractionPotential::MOLIERE, InteractionPotential::ZBL, InteractionPotential::LENZ_JENSEN] { From cc333157e4b1c308bb0474a76b8a8ddcd46aa5f2 Mon Sep 17 00:00:00 2001 From: Jon Drobny Date: Mon, 1 Mar 2021 10:45:50 -0600 Subject: [PATCH 09/14] Geometry rewrite completed. Next step will be testing Mesh0D to ensure that it works and reproduces equivalent results from Mesh2D. 3D targets are now possible (with a lot of dev work)! --- src/input.rs | 28 ++++++++++----- src/main.rs | 2 +- src/material.rs | 10 +++--- src/mesh.rs | 35 +++++++++++-------- src/tests.rs | 92 ++----------------------------------------------- 5 files changed, 46 insertions(+), 121 deletions(-) diff --git a/src/input.rs b/src/input.rs index b821750..86e6e63 100644 --- a/src/input.rs +++ b/src/input.rs @@ -1,16 +1,19 @@ use super::*; -pub trait InputFile { - type GeometryInput; +pub trait InputFile: GeometryInput { fn new(string: &str) -> Self; fn get_options(&self) -> &Options; fn get_material_parameters(&self) -> &material::MaterialParameters; fn get_particle_parameters(&self) -> &particle::ParticleParameters; - fn get_geometry_input(&self) -> &Self::GeometryInput; + fn get_geometry_input(&self) -> &::GeometryInput; +} + +pub trait GeometryInput { + type GeometryInput; } /// Rustbca's internal representation of an input file. -#[derive(Deserialize)] +#[derive(Deserialize, Clone)] pub struct Input2D { pub options: Options, pub material_parameters: material::MaterialParameters, @@ -18,8 +21,11 @@ pub struct Input2D { pub geometry_input: mesh::Mesh2DInput, } -impl InputFile for Input2D { +impl GeometryInput for Input2D { type GeometryInput = mesh::Mesh2DInput; +} + +impl InputFile for Input2D { fn new(string: &str) -> Input2D { toml::from_str(string).expect("Could not parse TOML file.") @@ -39,7 +45,7 @@ impl InputFile for Input2D { } } -#[derive(Deserialize)] +#[derive(Deserialize, Clone)] pub struct Input0D { pub options: Options, pub material_parameters: material::MaterialParameters, @@ -47,8 +53,11 @@ pub struct Input0D { pub geometry_input: mesh::Mesh0DInput, } -impl InputFile for Input0D { +impl GeometryInput for Input0D { type GeometryInput = mesh::Mesh0DInput; +} + +impl InputFile for Input0D { fn new(string: &str) -> Input0D { toml::from_str(string).expect("Could not parse TOML file.") @@ -131,7 +140,7 @@ pub struct Options { } pub fn input() -> (Vec, material::Material, Options, OutputUnits) -where ::InputFileFormat: Deserialize<'static> + 'static { +where ::InputFileFormat: Deserialize<'static> + 'static, ::GeometryInput: Clone { let args: Vec = env::args().collect(); @@ -154,9 +163,10 @@ where ::InputFileFormat: Deserialize<'static> + 'static { let input: ::InputFileFormat = InputFile::new(&input_toml); //Unpack toml information into structs - let material: material::Material = material::Material::::new(&input); let options = (*input.get_options()).clone(); let particle_parameters = (*input.get_particle_parameters()).clone(); + let material_parameters = (*input.get_material_parameters()).clone(); + let material: material::Material = material::Material::::new(material_parameters, input.get_geometry_input()); //Ensure nonsensical threads/chunks options crash on input assert!(options.num_threads > 0, "Input error: num_threads must be greater than zero."); diff --git a/src/main.rs b/src/main.rs index 428d002..d6886e8 100644 --- a/src/main.rs +++ b/src/main.rs @@ -62,7 +62,7 @@ pub mod structs; pub use crate::enums::*; pub use crate::consts::*; pub use crate::structs::*; -pub use crate::input::{Input2D, Input0D, Options, InputFile}; +pub use crate::input::{Input2D, Input0D, Options, InputFile, GeometryInput}; pub use crate::output::{OutputUnits}; pub use crate::mesh::{Geometry, GeometryElement}; diff --git a/src/material.rs b/src/material.rs index 694fdfd..cbc6d7d 100644 --- a/src/material.rs +++ b/src/material.rs @@ -1,7 +1,7 @@ use super::*; /// Holds material input parameters from [material_params]. -#[derive(Deserialize)] +#[derive(Deserialize, Clone)] pub struct MaterialParameters { pub energy_unit: String, pub mass_unit: String, @@ -26,11 +26,9 @@ pub struct Material { pub surface_binding_model: SurfaceBindingModel } -impl Material { +impl Material { - pub fn new(raw_input: &::InputFileFormat) -> Material { - - let material_parameters = raw_input.get_material_parameters(); + pub fn new(material_parameters: MaterialParameters, geometry_input: &<::InputFileFormat as GeometryInput>::GeometryInput) -> Material { let energy_unit: f64 = match material_parameters.energy_unit.as_str() { "EV" => EV, @@ -60,7 +58,7 @@ impl Material { Ec: material_parameters.Ec.iter().map(|&i| i*energy_unit).collect(), interaction_index: material_parameters.interaction_index.clone(), surface_binding_model: material_parameters.surface_binding_model, - geometry: Box::new(T::new(raw_input)), + geometry: Box::new(T::new(geometry_input)), } } diff --git a/src/mesh.rs b/src/mesh.rs index e585c0a..91346d6 100644 --- a/src/mesh.rs +++ b/src/mesh.rs @@ -3,14 +3,12 @@ use geo::algorithm::contains::Contains; use geo::{Polygon, LineString, Point, point, Closest}; use geo::algorithm::closest_point::ClosestPoint; - ///Trait for a Geometry object - all forms of geometry must implement these traits to be used -pub trait Geometry { - type GeometryInput; - type InputFileFormat: InputFile; +pub trait Geometry: GeometryInput { - fn new(input: &Self::InputFileFormat) -> Self; + type InputFileFormat: InputFile + Clone; + fn new(input: &<::InputFileFormat as GeometryInput>::GeometryInput) -> Self; fn get_densities(&self, x: f64, y: f64, z: f64) -> &Vec; fn get_ck(&self, x: f64, y: f64, z: f64) -> f64; fn get_total_density(&self, x: f64, y: f64, z: f64) -> f64; @@ -24,7 +22,7 @@ pub trait Geometry { fn closest_point(&self, x: f64, y: f64, z: f64) -> (f64, f64, f64); } -pub trait GeometryElement { +pub trait GeometryElement: Clone { fn contains(&self, x: f64, y: f64, z: f64) -> bool; fn distance_to(&self, x: f64, y: f64, z: f64) -> f64; fn get_densities(&self) -> &Vec; @@ -32,7 +30,7 @@ pub trait GeometryElement { fn get_electronic_stopping_correction_factor(&self) -> f64; } -#[derive(Deserialize)] +#[derive(Deserialize, Clone)] pub struct Mesh0DInput { pub length_unit: String, pub densities: Vec, @@ -40,6 +38,7 @@ pub struct Mesh0DInput { pub energy_barrier_thickness: f64, } +#[derive(Clone)] pub struct Mesh0D { pub densities: Vec, pub concentrations: Vec, @@ -47,14 +46,15 @@ pub struct Mesh0D { pub energy_barrier_thickness: f64, } +impl GeometryInput for Mesh0D { + type GeometryInput = Mesh0DInput; +} + impl Geometry for Mesh0D { - type GeometryInput = Mesh0DInput; type InputFileFormat = Input0D; - fn new(input: &Self::InputFileFormat) -> Mesh0D { - - let input = input.get_geometry_input(); + fn new(input: &::GeometryInput) -> Mesh0D { let length_unit: f64 = match input.length_unit.as_str() { "MICRON" => MICRON, @@ -123,7 +123,7 @@ impl Geometry for Mesh0D { } /// Object that contains raw mesh input data. -#[derive(Deserialize)] +#[derive(Deserialize, Clone)] pub struct Mesh2DInput { pub length_unit: String, pub triangles: Vec<(f64, f64, f64, f64, f64, f64)>, @@ -135,6 +135,7 @@ pub struct Mesh2DInput { } /// Triangular mesh for rustbca. +#[derive(Clone)] pub struct Mesh2D { mesh: Vec, pub boundary: Polygon, @@ -161,15 +162,17 @@ impl Mesh2D { } } +impl GeometryInput for Mesh2D { + type GeometryInput = Mesh2DInput; +} + impl Geometry for Mesh2D { type InputFileFormat = Input2D; - type GeometryInput = Mesh2DInput; /// Constructor for Mesh2D object from geometry_input. - fn new(input: &Self::InputFileFormat) -> Mesh2D { + fn new(geometry_input: &::GeometryInput) -> Mesh2D { - let geometry_input = input.get_geometry_input(); let triangles = geometry_input.triangles.clone(); let material_boundary_points = geometry_input.material_boundary_points.clone(); let simulation_boundary_points = geometry_input.simulation_boundary_points.clone(); @@ -329,6 +332,7 @@ impl Geometry for Mesh2D { } /// A mesh cell that contains a triangle and the local number densities and concentrations. +#[derive(Clone)] pub struct Cell2D { triangle: Triangle2D, pub densities: Vec, @@ -372,6 +376,7 @@ impl GeometryElement for Cell2D { } /// A triangle in 2D, with points (x1, y1), (x2, y2), (x3, y3), and the three line segments bewtween them. +#[derive(Clone)] pub struct Triangle2D { x1: f64, x2: f64, diff --git a/src/tests.rs b/src/tests.rs index ddf1fa5..fb54dc5 100644 --- a/src/tests.rs +++ b/src/tests.rs @@ -42,51 +42,7 @@ fn test_surface_binding_energy_barrier() { electronic_stopping_correction_factors: vec![0.0, 0.0], }; - let options = Options { - name: "test".to_string(), - track_trajectories: false, - track_recoils: true, - track_recoil_trajectories: false, - write_buffer_size: 8000, - weak_collision_order: 0, - suppress_deep_recoils: false, - high_energy_free_flight_paths: true, - electronic_stopping_mode: ElectronicStoppingMode::INTERPOLATED, - mean_free_path_model: MeanFreePathModel::LIQUID, - interaction_potential: vec![vec![InteractionPotential::MOLIERE]], - scattering_integral: vec![vec![ScatteringIntegral::MENDENHALL_WELLER]], - num_threads: 1, - num_chunks: 1, - use_hdf5: false, - root_finder: vec![vec![Rootfinder::NEWTON{max_iterations: 100, tolerance: 1E-3}]], - track_displacements: false, - track_energy_losses: false, - }; - - let particle_parameters = particle::ParticleParameters { - particle_input_filename: "test".to_string(), - length_unit: "ANGSTROM".to_string(), - energy_unit: "EV".to_string(), - mass_unit: "AMU".to_string(), - N: vec![], - m: vec![], - Z: vec![], - E: vec![], - Ec: vec![], - Es: vec![], - pos: vec![], - dir: vec![], - interaction_index: vec![], - }; - - let raw_input = Input2D { - options, - geometry_input, - material_parameters, - particle_parameters, - }; - - let material_1: material::Material = material::Material::::new(&raw_input); + let material_1: material::Material = material::Material::::new(material_parameters, &geometry_input); particle_1.pos.x = 500.*ANGSTROM; particle_1.pos.y = 0.; @@ -269,51 +225,7 @@ fn test_momentum_conservation() { energy_barrier_thickness: 0., }; - let options = Options { - name: "test".to_string(), - track_trajectories: false, - track_recoils: true, - track_recoil_trajectories: false, - write_buffer_size: 8000, - weak_collision_order: 0, - suppress_deep_recoils: false, - high_energy_free_flight_paths: true, - electronic_stopping_mode: ElectronicStoppingMode::INTERPOLATED, - mean_free_path_model: MeanFreePathModel::LIQUID, - interaction_potential: vec![vec![InteractionPotential::MOLIERE]], - scattering_integral: vec![vec![ScatteringIntegral::MENDENHALL_WELLER]], - num_threads: 1, - num_chunks: 1, - use_hdf5: false, - root_finder: vec![vec![Rootfinder::NEWTON{max_iterations: 100, tolerance: 1E-3}]], - track_displacements: false, - track_energy_losses: false, - }; - - let particle_parameters = particle::ParticleParameters { - particle_input_filename: "test".to_string(), - length_unit: "ANGSTROM".to_string(), - energy_unit: "EV".to_string(), - mass_unit: "AMU".to_string(), - N: vec![], - m: vec![], - Z: vec![], - E: vec![], - Ec: vec![], - Es: vec![], - pos: vec![], - dir: vec![], - interaction_index: vec![], - }; - - let raw_input = Input2D { - options, - geometry_input, - material_parameters, - particle_parameters, - }; - - let material_1: material::Material = material::Material::::new(&raw_input); + let material_1: material::Material = material::Material::::new(material_parameters, &geometry_input); for high_energy_free_flight_paths in vec![true, false] { for potential in vec![InteractionPotential::KR_C, InteractionPotential::MOLIERE, InteractionPotential::ZBL, InteractionPotential::LENZ_JENSEN] { From b71931aa574600bade32294428ed34076cf8a942 Mon Sep 17 00:00:00 2001 From: Jon Drobny Date: Mon, 1 Mar 2021 12:47:44 -0600 Subject: [PATCH 10/14] Apparently working version of 0D input file. --- examples/boron_nitride_0D.toml | 51 ++++++++++++++++++++++++++++++++++ src/bca.rs | 8 ++++++ src/main.rs | 3 ++ src/material.rs | 3 ++ src/mesh.rs | 9 +++++- 5 files changed, 73 insertions(+), 1 deletion(-) create mode 100644 examples/boron_nitride_0D.toml diff --git a/examples/boron_nitride_0D.toml b/examples/boron_nitride_0D.toml new file mode 100644 index 0000000..54fb8c9 --- /dev/null +++ b/examples/boron_nitride_0D.toml @@ -0,0 +1,51 @@ +[options] +name = "boron_dust_grain_" +track_trajectories = false +track_recoils = true +track_recoil_trajectories = false +write_buffer_size = 8000 +weak_collision_order = 0 +suppress_deep_recoils = false +high_energy_free_flight_paths = false +num_threads = 1 +num_chunks = 1 +use_hdf5 = false +electronic_stopping_mode = "LOW_ENERGY_NONLOCAL" +mean_free_path_model = "LIQUID" +interaction_potential = [["KR_C"]] +scattering_integral = [["MENDENHALL_WELLER"]] +root_finder = [[{"NEWTON"={max_iterations=100, tolerance=1E-3}}]] +track_displacements = false +track_energy_losses = false + +[material_parameters] +energy_unit = "EV" +mass_unit = "AMU" +Eb = [ 0.0, 0.0,] +Es = [ 5.76, 0.0,] +Ec = [ 1.0, 1.0,] +Z = [ 5, 7,] +m = [ 10.811, 14,] +interaction_index = [0, 0] +surface_binding_model = "AVERAGE" + +[particle_parameters] +length_unit = "MICRON" +energy_unit = "EV" +mass_unit = "AMU" +N = [ 10000 ] +m = [ 1.008 ] +Z = [ 1 ] +E = [ 1000.0 ] +Ec = [ 1.0 ] +Es = [ 10.0 ] +pos = [ [ 0.0, 0.0, 0.0,] ] +dir = [ [ 0.9999999999984769, 1.7453292519934434e-6, 0.0,] ] +interaction_index = [ 0 ] +particle_input_filename="" + +[geometry_input] +length_unit = "ANGSTROM" +energy_barrier_thickness = 1 +electronic_stopping_correction_factor = 1.0 +densities = [0.065, 0.065] diff --git a/src/bca.rs b/src/bca.rs index 96654c9..d1ed73f 100644 --- a/src/bca.rs +++ b/src/bca.rs @@ -76,6 +76,8 @@ pub fn single_ion_bca(particle: particle::Particle, material: &mate //Remove particle from top of vector as particle_1 let mut particle_1 = particles.pop().unwrap(); + //println!("Particle start Z: {} E: {} ({}, {}, {})", particle_1.Z, particle_1.E/Q, particle_1.pos.x/ANGSTROM, particle_1.pos.y/ANGSTROM, particle_1.pos.z/ANGSTROM); + //BCA loop while !particle_1.stopped & !particle_1.left { @@ -104,6 +106,8 @@ pub fn single_ion_bca(particle: particle::Particle, material: &mate particle_1.pos.x, particle_2.pos.x, &binary_collision_geometry)) .unwrap(); + //println!("{}", binary_collision_result); + //Only use 0th order collision for local electronic stopping if k == 0 { normalized_distance_of_closest_approach = binary_collision_result.normalized_distance_of_closest_approach; @@ -178,6 +182,9 @@ pub fn single_ion_bca(particle: particle::Particle, material: &mate bca::update_particle_energy(&mut particle_1, &material, distance_traveled, total_energy_loss, normalized_distance_of_closest_approach, strong_collision_Z, strong_collision_index, &options); + //println!("Particle finished collision loop Z: {} E: {} ({}, {}, {})", particle_1.Z, particle_1.E/Q, particle_1.pos.x/ANGSTROM, particle_1.pos.y/ANGSTROM, particle_1.pos.z/ANGSTROM); + + //println!("{} {} {}", energy_0/EV, energy_1/EV, (energy_1 - energy_0)/EV); //Check boundary conditions on leaving and stopping material::boundary_condition_planar(&mut particle_1, &material); @@ -185,6 +192,7 @@ pub fn single_ion_bca(particle: particle::Particle, material: &mate //Set particle index to topmost particle particle_index = particles.len(); } + //println!("Particle stopped or left Z: {} E: {} ({}, {}, {})", particle_1.Z, particle_1.E/Q, particle_1.pos.x/ANGSTROM, particle_1.pos.y/ANGSTROM, particle_1.pos.z/ANGSTROM); particle_output.push(particle_1); } particle_output diff --git a/src/main.rs b/src/main.rs index d6886e8..b6f4b22 100644 --- a/src/main.rs +++ b/src/main.rs @@ -68,7 +68,10 @@ pub use crate::mesh::{Geometry, GeometryElement}; fn main() { //Open and process input_file + #[cfg(not(feature = "mesh_0D"))] let (particle_input_array, material, options, output_units): (Vec, material::Material, Options, OutputUnits) = input::input(); + #[cfg(feature = "mesh_0D")] + let (particle_input_array, material, options, output_units): (Vec, material::Material, Options, OutputUnits) = input::input(); println!("Processing {} ions...", particle_input_array.len()); diff --git a/src/material.rs b/src/material.rs index cbc6d7d..38ecd2b 100644 --- a/src/material.rs +++ b/src/material.rs @@ -288,6 +288,7 @@ pub fn surface_binding_energy(particle_1: &mut particle::Particle, //Actual surface binding energies let Es = material.actual_surface_binding_energy(particle_1, x_old, y_old, z_old); + //println!("Actual Es: {}", Es); let Ec = particle_1.Ec; let inside_now = material.inside_energy_barrier(x, y, z); @@ -355,6 +356,8 @@ pub fn boundary_condition_planar(particle_1: &mut particle::Particl let E = particle_1.E; let Ec = particle_1.Ec; + //println!("({} {} {}) inside: {} inside energy: {} inside sim: {}", x/ANGSTROM, y/ANGSTROM, z/ANGSTROM, material.inside(x, y, z), material.inside_energy_barrier(x, y, z), material.inside_simulation_boundary(x, y, z)); + if !material.inside_simulation_boundary(x, y, z) { particle_1.left = true; particle_1.add_trajectory(); diff --git a/src/mesh.rs b/src/mesh.rs index 91346d6..bf4166a 100644 --- a/src/mesh.rs +++ b/src/mesh.rs @@ -70,13 +70,17 @@ impl Geometry for Mesh0D { }; let electronic_stopping_correction_factor = input.electronic_stopping_correction_factor; - let energy_barrier_thickness = input.energy_barrier_thickness/length_unit; + let energy_barrier_thickness = input.energy_barrier_thickness*length_unit; let densities: Vec = input.densities.iter().map(|element| element/(length_unit).powf(3.)).collect(); let total_density: f64 = densities.iter().sum(); let concentrations: Vec = densities.iter().map(|&density| density/total_density).collect::>(); + for density in densities.iter() { + println!("density: {}", density); + } + Mesh0D { densities, concentrations, @@ -106,9 +110,12 @@ impl Geometry for Mesh0D { fn inside(&self, x: f64, y: f64, z: f64) -> bool { x > 0.0 } + fn inside_simulation_boundary(&self, x: f64, y: f64, z: f64) -> bool { + //println!("x: {} energy_barrier_thickness: {}", x/ANGSTROM, self.energy_barrier_thickness/ANGSTROM); x > -2.*self.energy_barrier_thickness } + fn inside_energy_barrier(&self, x: f64, y: f64, z: f64) -> bool { x > -self.energy_barrier_thickness } From d7b4ae610a225eb6b880a864affda90dc9a848be Mon Sep 17 00:00:00 2001 From: Jon Drobny <37962344+drobnyjt@users.noreply.github.com> Date: Mon, 1 Mar 2021 12:54:31 -0600 Subject: [PATCH 11/14] Update rustbca_compile_check.yml added mesh_0D input file --- .github/workflows/rustbca_compile_check.yml | 1 + 1 file changed, 1 insertion(+) diff --git a/.github/workflows/rustbca_compile_check.yml b/.github/workflows/rustbca_compile_check.yml index 17ed2de..56a8169 100644 --- a/.github/workflows/rustbca_compile_check.yml +++ b/.github/workflows/rustbca_compile_check.yml @@ -43,3 +43,4 @@ jobs: run: | cargo run --release examples/boron_nitride.toml cargo run --release --features distributions examples/layered_geometry.toml + cargo run --release --features mesh_0D examples/boron_nitride_0D.toml From 6d3d04d735ad10b2d27aaeb9062dcbf46bbab51f Mon Sep 17 00:00:00 2001 From: Jon Drobny Date: Mon, 1 Mar 2021 13:06:11 -0600 Subject: [PATCH 12/14] Forgot to add Cargo that includes the mesh_0D feature. --- Cargo.toml | 1 + 1 file changed, 1 insertion(+) diff --git a/Cargo.toml b/Cargo.toml index c004727..21a6283 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -39,3 +39,4 @@ cpr_rootfinder_netlib = ["rcpr", "netlib-src"] cpr_rootfinder_intel_mkl = ["rcpr", "intel-mkl-src"] distributions = ["ndarray"] no_list_output = [] +mesh_0D = [] From bac84bd2c11facde96d76522088a552268371b85 Mon Sep 17 00:00:00 2001 From: Jon Drobny Date: Mon, 1 Mar 2021 16:22:24 -0600 Subject: [PATCH 13/14] Mesh0D reproduces sputtering yields when compared to Mesh2D. Added some unit tests. --- examples/boron_nitride_0D.toml | 1 - examples/titanium_dioxide_0D.toml | 65 +++++++++++++ src/input.rs | 2 +- src/material.rs | 4 +- src/mesh.rs | 7 +- src/tests.rs | 147 +++++++++++++++++++++++++----- 6 files changed, 197 insertions(+), 29 deletions(-) create mode 100644 examples/titanium_dioxide_0D.toml diff --git a/examples/boron_nitride_0D.toml b/examples/boron_nitride_0D.toml index 54fb8c9..4aad3ca 100644 --- a/examples/boron_nitride_0D.toml +++ b/examples/boron_nitride_0D.toml @@ -46,6 +46,5 @@ particle_input_filename="" [geometry_input] length_unit = "ANGSTROM" -energy_barrier_thickness = 1 electronic_stopping_correction_factor = 1.0 densities = [0.065, 0.065] diff --git a/examples/titanium_dioxide_0D.toml b/examples/titanium_dioxide_0D.toml new file mode 100644 index 0000000..ff1e7f3 --- /dev/null +++ b/examples/titanium_dioxide_0D.toml @@ -0,0 +1,65 @@ +[options] +name = "2000.0eV_0.0001deg_He_TiO2" +track_trajectories = false +track_recoils = true +track_recoil_trajectories = false +write_buffer_size = 8000 +weak_collision_order = 3 +suppress_deep_recoils = false +high_energy_free_flight_paths = false +electronic_stopping_mode = "LOW_ENERGY_NONLOCAL" +mean_free_path_model = "LIQUID" +interaction_potential = [["KR_C"]] +scattering_integral = [["MENDENHALL_WELLER"]] +use_hdf5 = false +root_finder = [[{"NEWTON"={max_iterations=100, tolerance=1E-6}}]] +num_threads = 4 +num_chunks = 10 +track_displacements = false +track_energy_losses = false +energy_min = 0.0 +energy_max = 2000.0 +energy_num = 20 +angle_min = 0.0 +angle_max = 90.0 +angle_num = 20 +x_min = 0.0 +y_min = -0.1 +z_min = -0.1 +x_max = 0.1 +y_max = 0.1 +z_max = 0.1 +x_num = 20 +y_num = 20 +z_num = 20 + +[particle_parameters] +length_unit = "1e-6" +energy_unit = "EV" +mass_unit = "AMU" +N = [ 100000,] +m = [ 4.002602,] +Z = [ 2,] +E = [ 2000.0,] +Ec = [ 1.0,] +Es = [ 0.0,] +interaction_index = [0] +pos = [ [ 0.0, 0.0, 0.0,],] +dir = [ [ 0.5, 0.5, 0.0,],] +particle_input_filename = "" + +[geometry_input] +length_unit = "MICRON" +densities = [ 3.0e+10, 6.0e+10, 0.0, 0.0 ] +electronic_stopping_correction_factor = 1.0 + +[material_parameters] +energy_unit = "EV" +mass_unit = "AMU" +Eb = [ 3.0, 2.58, 3.0, 0.0,] +Es = [ 4.84, 2.58, 3.39, 4.72,] +Ec = [ 3.5, 2.0, 3.0, 1.5,] +Z = [ 22, 8, 13, 14,] +m = [ 47.867, 15.9994, 26.98, 28.08553,] +interaction_index = [0, 0, 0, 0] +surface_binding_model = "TARGET" diff --git a/src/input.rs b/src/input.rs index 86e6e63..7315d4d 100644 --- a/src/input.rs +++ b/src/input.rs @@ -166,7 +166,7 @@ where ::InputFileFormat: Deserialize<'static> + 'static, = material::Material::::new(material_parameters, input.get_geometry_input()); + let material: material::Material = material::Material::::new(&material_parameters, input.get_geometry_input()); //Ensure nonsensical threads/chunks options crash on input assert!(options.num_threads > 0, "Input error: num_threads must be greater than zero."); diff --git a/src/material.rs b/src/material.rs index 38ecd2b..4c110f0 100644 --- a/src/material.rs +++ b/src/material.rs @@ -28,7 +28,7 @@ pub struct Material { } impl Material { - pub fn new(material_parameters: MaterialParameters, geometry_input: &<::InputFileFormat as GeometryInput>::GeometryInput) -> Material { + pub fn new(material_parameters: &MaterialParameters, geometry_input: &<::InputFileFormat as GeometryInput>::GeometryInput) -> Material { let energy_unit: f64 = match material_parameters.energy_unit.as_str() { "EV" => EV, @@ -356,8 +356,6 @@ pub fn boundary_condition_planar(particle_1: &mut particle::Particl let E = particle_1.E; let Ec = particle_1.Ec; - //println!("({} {} {}) inside: {} inside energy: {} inside sim: {}", x/ANGSTROM, y/ANGSTROM, z/ANGSTROM, material.inside(x, y, z), material.inside_energy_barrier(x, y, z), material.inside_simulation_boundary(x, y, z)); - if !material.inside_simulation_boundary(x, y, z) { particle_1.left = true; particle_1.add_trajectory(); diff --git a/src/mesh.rs b/src/mesh.rs index bf4166a..1fa9619 100644 --- a/src/mesh.rs +++ b/src/mesh.rs @@ -35,7 +35,6 @@ pub struct Mesh0DInput { pub length_unit: String, pub densities: Vec, pub electronic_stopping_correction_factor: f64, - pub energy_barrier_thickness: f64, } #[derive(Clone)] @@ -70,11 +69,13 @@ impl Geometry for Mesh0D { }; let electronic_stopping_correction_factor = input.electronic_stopping_correction_factor; - let energy_barrier_thickness = input.energy_barrier_thickness*length_unit; let densities: Vec = input.densities.iter().map(|element| element/(length_unit).powf(3.)).collect(); let total_density: f64 = densities.iter().sum(); + + let energy_barrier_thickness = total_density.powf(-1./3.)/SQRTPI*2.; + let concentrations: Vec = densities.iter().map(|&density| density/total_density).collect::>(); for density in densities.iter() { @@ -113,7 +114,7 @@ impl Geometry for Mesh0D { fn inside_simulation_boundary(&self, x: f64, y: f64, z: f64) -> bool { //println!("x: {} energy_barrier_thickness: {}", x/ANGSTROM, self.energy_barrier_thickness/ANGSTROM); - x > -2.*self.energy_barrier_thickness + x > -10.*self.energy_barrier_thickness } fn inside_energy_barrier(&self, x: f64, y: f64, z: f64) -> bool { diff --git a/src/tests.rs b/src/tests.rs index fb54dc5..d6eaebc 100644 --- a/src/tests.rs +++ b/src/tests.rs @@ -3,6 +3,72 @@ use super::*; #[cfg(test)] use float_cmp::*; +#[test] +fn test_geometry() { + let mass = 1.; + let Z = 1.; + let E = 10.*EV; + let Ec = 1.*EV; + let Es = 5.76*EV; + let x = 0.; + let y = 0.; + let z = 0.; + let cosx = 1./(2.0_f64).sqrt(); + let cosy = 1./(2.0_f64).sqrt(); + let cosz = 0.; + let mut particle_1 = particle::Particle::new(mass, Z, E, Ec, Es, x, y, z, cosx, cosy, cosz, false, false, 0); + + let material_parameters = material::MaterialParameters{ + energy_unit: "EV".to_string(), + mass_unit: "AMU".to_string(), + Eb: vec![0.0, 0.0], + Es: vec![2.0, 4.0], + Ec: vec![1.0, 1.0], + Z: vec![29., 1.], + m: vec![63.54, 1.0008], + interaction_index: vec![0, 0], + surface_binding_model: SurfaceBindingModel::TARGET + }; + + let thickness: f64 = 1000.; + let depth: f64 = 1000.; + + let geometry_input_2D = mesh::Mesh2DInput { + length_unit: "ANGSTROM".to_string(), + triangles: vec![(0., depth, 0., thickness/2., thickness/2., -thickness/2.), (depth, depth, 0., thickness/2., -thickness/2., -thickness/2.)], + densities: vec![vec![0.03, 0.03], vec![0.03, 0.03]], + material_boundary_points: vec![(0., thickness/2.), (depth, thickness/2.), (depth, -thickness/2.), (0., -thickness/2.), (0., thickness/2.)], + simulation_boundary_points: vec![(0., 1.1*thickness/2.), (depth, 1.1*thickness/2.), (depth, -1.1*thickness/2.), (0., -1.1*thickness/2.), (0., 1.1*thickness/2.)], + energy_barrier_thickness: 10., + electronic_stopping_correction_factors: vec![1.0, 1.0], + }; + + let geometry_input_0D = mesh::Mesh0DInput { + length_unit: "ANGSTROM".to_string(), + densities: vec![0.03, 0.03], + electronic_stopping_correction_factor: 1.0 + }; + + let material_2D: material::Material = material::Material::::new(&material_parameters, &geometry_input_2D); + let mut material_0D: material::Material = material::Material::::new(&material_parameters, &geometry_input_0D); + material_0D.geometry.energy_barrier_thickness = 10.*ANGSTROM; + + particle_1.pos.x = 500.*ANGSTROM; + particle_1.pos.y = 0.; + + particle_1.pos_old.x = -500.*ANGSTROM; + particle_1.pos_old.y = 0.; + + assert_eq!(material_2D.geometry.get_ck(0., 0., 0.), material_0D.geometry.get_ck(0., 0., 0.)); + assert_eq!(material_2D.geometry.get_energy_barrier_thickness(), material_0D.geometry.get_energy_barrier_thickness()); + assert_eq!(material_2D.geometry.get_densities(0., 0., 0.), material_0D.geometry.get_densities(0., 0., 0.)); + assert_eq!(material_2D.geometry.get_total_density(0., 0., 0.), material_0D.geometry.get_total_density(0., 0., 0.)); + assert_eq!(material_2D.geometry.get_concentrations(0., 0., 0.), material_0D.geometry.get_concentrations(0., 0., 0.)); + assert_eq!(material_2D.geometry.closest_point(-10., 0., 5.), material_0D.geometry.closest_point(-10., 0., 5.)); + assert_eq!(material_2D.geometry.get_densities_nearest_to(-10., 0., 5.), material_0D.geometry.get_densities_nearest_to(-10., 0., 5.)); + assert_eq!(material_2D.geometry.get_ck_nearest_to(-10., 0., 5.), material_0D.geometry.get_ck_nearest_to(-10., 0., 5.)); +} + #[test] fn test_surface_binding_energy_barrier() { let mass = 1.; @@ -32,17 +98,26 @@ fn test_surface_binding_energy_barrier() { let thickness: f64 = 1000.; let depth: f64 = 1000.; - let geometry_input = mesh::Mesh2DInput { + + let geometry_input_2D = mesh::Mesh2DInput { length_unit: "ANGSTROM".to_string(), triangles: vec![(0., depth, 0., thickness/2., thickness/2., -thickness/2.), (depth, depth, 0., thickness/2., -thickness/2., -thickness/2.)], densities: vec![vec![0.03, 0.03], vec![0.03, 0.03]], material_boundary_points: vec![(0., thickness/2.), (depth, thickness/2.), (depth, -thickness/2.), (0., -thickness/2.), (0., thickness/2.)], simulation_boundary_points: vec![(0., 1.1*thickness/2.), (depth, 1.1*thickness/2.), (depth, -1.1*thickness/2.), (0., -1.1*thickness/2.), (0., 1.1*thickness/2.)], energy_barrier_thickness: 10., - electronic_stopping_correction_factors: vec![0.0, 0.0], + electronic_stopping_correction_factors: vec![1.0, 1.0], + }; + + let geometry_input_0D = mesh::Mesh0DInput { + length_unit: "ANGSTROM".to_string(), + densities: vec![0.03, 0.03], + electronic_stopping_correction_factor: 1.0 }; - let material_1: material::Material = material::Material::::new(material_parameters, &geometry_input); + let material_2D: material::Material = material::Material::::new(&material_parameters, &geometry_input_2D); + let mut material_0D: material::Material = material::Material::::new(&material_parameters, &geometry_input_0D); + material_0D.geometry.energy_barrier_thickness = 10.*ANGSTROM; particle_1.pos.x = 500.*ANGSTROM; particle_1.pos.y = 0.; @@ -50,35 +125,65 @@ fn test_surface_binding_energy_barrier() { particle_1.pos_old.x = -500.*ANGSTROM; particle_1.pos_old.y = 0.; - let inside = material_1.inside(particle_1.pos.x, particle_1.pos.y, particle_1.pos.z); - let inside_old = material_1.inside(particle_1.pos_old.x, particle_1.pos_old.y, particle_1.pos.z); + let inside_0D = material_0D.inside(particle_1.pos.x, particle_1.pos.y, particle_1.pos.z); + let inside_old_0D = material_0D.inside(particle_1.pos_old.x, particle_1.pos_old.y, particle_1.pos_old.z); + + let inside_2D = material_2D.inside(particle_1.pos.x, particle_1.pos.y, particle_1.pos.z); + let inside_old_2D = material_2D.inside(particle_1.pos_old.x, particle_1.pos_old.y, particle_1.pos_old.z); //println!("{} {}", inside, inside_old); - assert!(inside); - assert!(!inside_old); + assert!(inside_0D); + assert!(!inside_old_0D); + + assert!(inside_2D); + assert!(!inside_old_2D); //Test concentration-dependent surface binding energy - let surface_binding_energy = material_1.actual_surface_binding_energy(&particle_1, particle_1.pos.x, particle_1.pos.y, particle_1.pos.z); - assert!(approx_eq!(f64, surface_binding_energy/EV, (2. + 4.)/2., epsilon=1E-12)); - //println!("sbv: {}", surface_binding_energy/EV); + let surface_binding_energy_2D = material_2D.actual_surface_binding_energy(&particle_1, particle_1.pos.x, particle_1.pos.y, particle_1.pos.z); + let surface_binding_energy_0D = material_0D.actual_surface_binding_energy(&particle_1, particle_1.pos.x, particle_1.pos.y, particle_1.pos.z); + + assert!(approx_eq!(f64, surface_binding_energy_0D/EV, (2. + 4.)/2., epsilon=1E-24)); + assert!(approx_eq!(f64, surface_binding_energy_2D/EV, (2. + 4.)/2., epsilon=1E-24)); + + println!("sbv 0D: {}", surface_binding_energy_0D/EV); + println!("sbv 2D: {}", surface_binding_energy_2D/EV); //Test leftmost boundary - assert!(material_1.inside_energy_barrier(500.*ANGSTROM, 0., 0.)); - assert!(material_1.inside_energy_barrier(-5.*ANGSTROM, 0., 0.)); - assert!(!material_1.inside_energy_barrier(-15.*ANGSTROM, 0., 0.)); + assert!(material_2D.inside_energy_barrier(500.*ANGSTROM, 0., 0.)); + assert!(material_2D.inside_energy_barrier(-5.*ANGSTROM, 0., 0.)); + assert!(!material_2D.inside_energy_barrier(-15.*ANGSTROM, 0., 0.)); + + assert!(material_0D.inside_energy_barrier(500.*ANGSTROM, 0., 0.)); + assert!(material_0D.inside_energy_barrier(-5.*ANGSTROM, 0., 0.)); + assert!(!material_0D.inside_energy_barrier(-15.*ANGSTROM, 0., 0.)); + + assert!(material_0D.inside_energy_barrier(500.*ANGSTROM, 0., 0.)); + assert!(!material_0D.inside_energy_barrier(-500.*ANGSTROM, 0., 0.)); + assert!(material_0D.inside_energy_barrier(-9.*ANGSTROM, 0., 0.)); + assert!(!material_0D.inside_energy_barrier(-11.*ANGSTROM, 0., 0.)); //Test top boundary - assert!(material_1.inside_energy_barrier(500.*ANGSTROM, 0., 0.)); - assert!(material_1.inside_energy_barrier(500.*ANGSTROM, 505.*ANGSTROM, 0.)); - assert!(!material_1.inside_energy_barrier(500.*ANGSTROM, 515.*ANGSTROM, 0.)); + assert!(material_2D.inside_energy_barrier(500.*ANGSTROM, 0., 0.)); + assert!(material_2D.inside_energy_barrier(500.*ANGSTROM, 505.*ANGSTROM, 0.)); + assert!(!material_2D.inside_energy_barrier(500.*ANGSTROM, 515.*ANGSTROM, 0.)); + + assert!(material_0D.inside_energy_barrier(500.*ANGSTROM, 0., 0.)); + assert!(material_0D.inside_energy_barrier(500.*ANGSTROM, 505.*ANGSTROM, 0.)); + assert!(material_0D.inside_energy_barrier(500.*ANGSTROM, 515.*ANGSTROM, 0.)); //Test bottom boundary - assert!(material_1.inside_energy_barrier(500.*ANGSTROM, -505.*ANGSTROM, 0.)); - assert!(!material_1.inside_energy_barrier(500.*ANGSTROM, -515.*ANGSTROM, 0.)); + assert!(material_2D.inside_energy_barrier(500.*ANGSTROM, -505.*ANGSTROM, 0.)); + assert!(!material_2D.inside_energy_barrier(500.*ANGSTROM, -515.*ANGSTROM, 0.)); + + assert!(material_0D.inside_energy_barrier(500.*ANGSTROM, -505.*ANGSTROM, 0.)); + assert!(material_0D.inside_energy_barrier(500.*ANGSTROM, -515.*ANGSTROM, 0.)); //Test rightmost boundary - assert!(material_1.inside_energy_barrier(1005.*ANGSTROM, 0., 0.)); - assert!(!material_1.inside_energy_barrier(1015.*ANGSTROM, 0., 0.)); + assert!(material_2D.inside_energy_barrier(1005.*ANGSTROM, 0., 0.)); + assert!(!material_2D.inside_energy_barrier(1015.*ANGSTROM, 0., 0.)); + + assert!(material_0D.inside_energy_barrier(1005.*ANGSTROM, 0., 0.)); + assert!(material_0D.inside_energy_barrier(1015.*ANGSTROM, 0., 0.)); } #[test] @@ -225,7 +330,7 @@ fn test_momentum_conservation() { energy_barrier_thickness: 0., }; - let material_1: material::Material = material::Material::::new(material_parameters, &geometry_input); + let material_1: material::Material = material::Material::::new(&material_parameters, &geometry_input); for high_energy_free_flight_paths in vec![true, false] { for potential in vec![InteractionPotential::KR_C, InteractionPotential::MOLIERE, InteractionPotential::ZBL, InteractionPotential::LENZ_JENSEN] { From 78fb8fb1c2b3c76f049234ff72d50a1dbe93abde Mon Sep 17 00:00:00 2001 From: Jon Drobny Date: Tue, 2 Mar 2021 15:22:29 -0600 Subject: [PATCH 14/14] Rewrite of input. Now, if you want to use 0D, you don't have to build with the 0D feature, you just give the code an additional flag. --- .github/workflows/rustbca_compile_check.yml | 5 +- Cargo.toml | 6 +- src/enums.rs | 11 ++ src/input.rs | 11 +- src/main.rs | 167 +++++++++++--------- src/material.rs | 1 - src/mesh.rs | 5 +- 7 files changed, 118 insertions(+), 88 deletions(-) diff --git a/.github/workflows/rustbca_compile_check.yml b/.github/workflows/rustbca_compile_check.yml index 56a8169..427930b 100644 --- a/.github/workflows/rustbca_compile_check.yml +++ b/.github/workflows/rustbca_compile_check.yml @@ -26,7 +26,7 @@ jobs: sudo apt-get install python3-pip - name: Install Python libraries run: | - pip3 + pip3 python3 -m pip install numpy shapely scipy - name: Install python TOML library from source run: | @@ -42,5 +42,6 @@ jobs: - name: Run Examples run: | cargo run --release examples/boron_nitride.toml + cargo run --release 0D examples/boron_nitride_0D.toml + cargo run --release 0D examples/titanium_dioxide_0D.toml cargo run --release --features distributions examples/layered_geometry.toml - cargo run --release --features mesh_0D examples/boron_nitride_0D.toml diff --git a/Cargo.toml b/Cargo.toml index 21a6283..53ee5ff 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -9,12 +9,13 @@ name="rustBCA" [dependencies] rand = "0.8.3" -geo = "0.17.1" + toml = "0.5.8" anyhow = "1.0.38" itertools = "0.10.0" rayon = "1.5.0" -indicatif = {version="0.15.0", features=["rayon"]} +geo = {version = "0.17.1", optional = false} +indicatif = {version = "0.15.0", features=["rayon"]} serde = { version = "1.0.123", features = ["derive"] } hdf5 = {version = "0.7.1", optional = true} openblas-src = {version = "0.9", optional = true} @@ -40,3 +41,4 @@ cpr_rootfinder_intel_mkl = ["rcpr", "intel-mkl-src"] distributions = ["ndarray"] no_list_output = [] mesh_0D = [] +mesh_2D = [] diff --git a/src/enums.rs b/src/enums.rs index a1dca03..0da78e5 100644 --- a/src/enums.rs +++ b/src/enums.rs @@ -1,5 +1,16 @@ use super::*; +pub enum MaterialType { + MESH0D(material::Material), + MESH2D(material::Material) +} + +#[derive(Deserialize)] +pub enum GeometryType { + MESH0D, + MESH2D, +} + /// Mode of electronic stopping to use. #[derive(Deserialize, PartialEq, Clone, Copy)] pub enum ElectronicStoppingMode { diff --git a/src/input.rs b/src/input.rs index 7315d4d..4188b1c 100644 --- a/src/input.rs +++ b/src/input.rs @@ -28,7 +28,7 @@ impl GeometryInput for Input2D { impl InputFile for Input2D { fn new(string: &str) -> Input2D { - toml::from_str(string).expect("Could not parse TOML file.") + toml::from_str(string).unwrap() } fn get_options(&self) -> &Options{ @@ -144,10 +144,11 @@ where ::InputFileFormat: Deserialize<'static> + 'static, = env::args().collect(); - let input_file = match args.len() { - 1 => "input.toml".to_string(), - 2 => args[1].clone(), - _ => panic!("Too many command line arguments. RustBCA accepts 0 (use 'input.toml') or 1 (input file name).") + let (input_file, geometry_type) = match args.len() { + 1 => ("input.toml".to_string(), GeometryType::MESH2D), + 2 => (args[1].clone(), GeometryType::MESH2D), + 3 => (args[2].clone(), match args[1].as_str() { "0D" => GeometryType::MESH0D, "2D" => GeometryType::MESH2D, _ => panic!("Unimplemented geometry {}.", args[1].clone()) }), + _ => panic!("Too many command line arguments. RustBCA accepts 0 (use 'input.toml') 1 () or 2 ( )"), }; //Read input file, convert to string, and open with toml diff --git a/src/main.rs b/src/main.rs index b6f4b22..c34d906 100644 --- a/src/main.rs +++ b/src/main.rs @@ -64,87 +64,106 @@ pub use crate::consts::*; pub use crate::structs::*; pub use crate::input::{Input2D, Input0D, Options, InputFile, GeometryInput}; pub use crate::output::{OutputUnits}; -pub use crate::mesh::{Geometry, GeometryElement}; +pub use crate::mesh::{Geometry, GeometryElement, Mesh0D, Mesh2D}; -fn main() { - //Open and process input_file - #[cfg(not(feature = "mesh_0D"))] - let (particle_input_array, material, options, output_units): (Vec, material::Material, Options, OutputUnits) = input::input(); - #[cfg(feature = "mesh_0D")] - let (particle_input_array, material, options, output_units): (Vec, material::Material, Options, OutputUnits) = input::input(); - - println!("Processing {} ions...", particle_input_array.len()); - - 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."); - - #[cfg(not(feature = "no_list_output"))] - let mut output_list_streams = output::open_output_lists(&options); - - let mut summary = output::SummaryPerSpecies::new(&options); - - #[cfg(feature = "distributions")] - let mut distributions = output::Distributions::new(&options); - - //Initialize threads with rayon - println!("Initializing with {} threads...", options.num_threads); - if options.num_threads > 1 {let pool = rayon::ThreadPoolBuilder::new().num_threads(options.num_threads).build_global().unwrap();}; - - //Create and configure progress bar - let bar: ProgressBar = ProgressBar::new(total_count); - bar.set_style(ProgressStyle::default_bar() - .template("[{elapsed_precise}][{bar:40.cyan/blue}][{eta_precise}] {percent}%") - .progress_chars("#>-")); - - //Main loop - for (chunk_index, particle_input_chunk) in particle_input_array.chunks((total_count/options.num_chunks) as usize).enumerate() { - - let mut finished_particles: Vec = Vec::new(); - - if options.num_threads > 1 { - // BCA loop is implemented as parallelized extension of a per-chunk initially empty - // finished particle array via map from particle -> finished particles via BCA - finished_particles.par_extend( - particle_input_chunk.into_par_iter() - .map(|particle_input| { - bar.tick(); - bar.inc(1); - bca::single_ion_bca(particle::Particle::from_input(*particle_input, &options), &material, &options) - }).flatten() - ); - } else { - finished_particles.extend( - particle_input_chunk.iter() - .map(|particle_input| { - bar.tick(); - bar.inc(1); - bca::single_ion_bca(particle::Particle::from_input(*particle_input, &options), &material, &options) - }).flatten() - ); - } +fn physics_loop(particle_input_array: Vec, material: material::Material, options: Options, output_units: OutputUnits) { - // Process this chunk of finished particles for output - for particle in finished_particles { + println!("Processing {} ions...", particle_input_array.len()); - summary.update(&particle); - - #[cfg(feature = "distributions")] - distributions.update(&particle, &output_units); + 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."); + #[cfg(not(feature = "no_list_output"))] + let mut output_list_streams = output::open_output_lists(&options); + + let mut summary = output::SummaryPerSpecies::new(&options); + + #[cfg(feature = "distributions")] + let mut distributions = output::Distributions::new(&options); + + //Initialize threads with rayon + println!("Initializing with {} threads...", options.num_threads); + if options.num_threads > 1 {let pool = rayon::ThreadPoolBuilder::new().num_threads(options.num_threads).build_global().unwrap();}; + + //Create and configure progress bar + let bar: ProgressBar = ProgressBar::new(total_count); + bar.set_style(ProgressStyle::default_bar() + .template("[{elapsed_precise}][{bar:40.cyan/blue}][{eta_precise}] {percent}%") + .progress_chars("#>-")); + + //Main loop + for (chunk_index, particle_input_chunk) in particle_input_array.chunks((total_count/options.num_chunks) as usize).enumerate() { + + let mut finished_particles: Vec = Vec::new(); + + if options.num_threads > 1 { + // BCA loop is implemented as parallelized extension of a per-chunk initially empty + // finished particle array via map from particle -> finished particles via BCA + finished_particles.par_extend( + particle_input_chunk.into_par_iter() + .map(|particle_input| { + bar.tick(); + bar.inc(1); + bca::single_ion_bca(particle::Particle::from_input(*particle_input, &options), &material, &options) + }).flatten() + ); + } else { + finished_particles.extend( + particle_input_chunk.iter() + .map(|particle_input| { + bar.tick(); + bar.inc(1); + bca::single_ion_bca(particle::Particle::from_input(*particle_input, &options), &material, &options) + }).flatten() + ); + } + + // Process this chunk of finished particles for output + for particle in finished_particles { + + summary.update(&particle); + + #[cfg(feature = "distributions")] + distributions.update(&particle, &output_units); + + #[cfg(not(feature = "no_list_output"))] + output::output_lists(&mut output_list_streams, particle, &options, &output_units); + + } + //Flush all file streams before dropping to ensure all data is written #[cfg(not(feature = "no_list_output"))] - output::output_lists(&mut output_list_streams, particle, &options, &output_units); - + output::output_list_flush(&mut output_list_streams); } - //Flush all file streams before dropping to ensure all data is written - #[cfg(not(feature = "no_list_output"))] - output::output_list_flush(&mut output_list_streams); - } - summary.print(&options, &output_units); + summary.print(&options, &output_units); - //Write distributions to file - #[cfg(feature = "distributions")] - distributions.print(&options); + //Write distributions to file + #[cfg(feature = "distributions")] + distributions.print(&options); + + println!("Finished!"); +} + +fn main() { + + let args: Vec = env::args().collect(); + + let (input_file, geometry_type) = match args.len() { + 1 => ("input.toml".to_string(), GeometryType::MESH2D), + 2 => (args[1].clone(), GeometryType::MESH2D), + 3 => (args[2].clone(), match args[1].as_str() { "0D" => GeometryType::MESH0D, "2D" => GeometryType::MESH2D, _ => panic!("Unimplemented geometry {}.", args[1].clone()) }), + _ => 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::(); + physics_loop::(particle_input_array, material, options, output_units); + }, + GeometryType::MESH2D => { + let (particle_input_array, material, options, output_units) = input::input::(); + physics_loop::(particle_input_array, material, options, output_units); + }, + } - println!("Finished!"); } diff --git a/src/material.rs b/src/material.rs index 4c110f0..174f3e0 100644 --- a/src/material.rs +++ b/src/material.rs @@ -24,7 +24,6 @@ pub struct Material { pub interaction_index: Vec, pub geometry: Box, pub surface_binding_model: SurfaceBindingModel - } impl Material { diff --git a/src/mesh.rs b/src/mesh.rs index 1fa9619..6985696 100644 --- a/src/mesh.rs +++ b/src/mesh.rs @@ -1,4 +1,5 @@ use super::*; + use geo::algorithm::contains::Contains; use geo::{Polygon, LineString, Point, point, Closest}; use geo::algorithm::closest_point::ClosestPoint; @@ -78,10 +79,6 @@ impl Geometry for Mesh0D { let concentrations: Vec = densities.iter().map(|&density| density/total_density).collect::>(); - for density in densities.iter() { - println!("density: {}", density); - } - Mesh0D { densities, concentrations,