From 34f9318cce690003029b08cf24a766320a01d439 Mon Sep 17 00:00:00 2001 From: Sam Morrell Date: Wed, 21 Feb 2024 16:45:47 +0000 Subject: [PATCH 01/16] Initial implementation of boundary condition struct. - Added structs for Boundary, BoundaryCondition, BoundaryDirection. - Added handling of grid boundary hitting event. - Added hard-coded default behaviour to all engines. --- src/bin/mcrt.rs | 13 +- src/geom/domain/boundary.rs | 501 +++++++++++++++++++++++++ src/geom/domain/mod.rs | 5 +- src/sim/engine/engines/fluorescence.rs | 7 +- src/sim/engine/engines/photo.rs | 7 +- src/sim/engine/engines/raman.rs | 7 +- src/sim/engine/engines/standard.rs | 9 +- src/sim/event.rs | 28 +- src/sim/input.rs | 7 +- 9 files changed, 569 insertions(+), 15 deletions(-) create mode 100644 src/geom/domain/boundary.rs diff --git a/src/bin/mcrt.rs b/src/bin/mcrt.rs index c0a31b6..8cb7394 100644 --- a/src/bin/mcrt.rs +++ b/src/bin/mcrt.rs @@ -10,8 +10,9 @@ use aetherus::{ args, data::Histogram, fs::{File, Load, Save}, - geom::{Grid, Tree}, + geom::{Grid, Tree, Boundary, BoundaryCondition}, img::{Colour, Image}, + math::{Point3}, ord::{Build, Link, Register, Set, cartesian::{X, Y}}, report, sim::{ @@ -84,6 +85,14 @@ fn main() { .expect("Failed to link attribute to surfaces."); report!(surfs, "surfaces"); + /* + * Create a boundary for the simulation with boundary coneditions. + * For now we hard-code this to kill, but we can link this to configuration soon. + * TODO: We probably want to implement the MPI adjacent rank transfer here too. + */ + let bound = Boundary::new_kill(grid.boundary().clone()); + report!(bound, "boundary conditions"); + sub_section(term_width, "Growing"); let tree = Tree::new(¶ms.tree, &surfs); report!(tree, "hit-scan tree"); @@ -95,7 +104,7 @@ fn main() { .fold(base_output.clone(), |mut output, (light_idx, (light_id, light))| { section(term_width, &format!("Running for light {} ({} / {})", light_id, light_idx + 1, nlights)); report!(light, light_id); - let input = Input::new(&spec_reg, &mats, &attrs, light, &tree, &grid, &sett); + let input = Input::new(&spec_reg, &mats, &attrs, light, &tree, &grid, &sett, &bound); let data = run::multi_thread(&engine, input, &base_output).expect("Failed to run MCRT."); diff --git a/src/geom/domain/boundary.rs b/src/geom/domain/boundary.rs new file mode 100644 index 0000000..328aba3 --- /dev/null +++ b/src/geom/domain/boundary.rs @@ -0,0 +1,501 @@ +use crate::{ + access, clone, fmt_report, geom::{Cube, Hit, Ray, Side, Trace}, math::{Dir3, Vec3}, phys::{Photon, Reflectance}, sim::Attribute +}; +use rand::rngs::ThreadRng; +use std::fmt::{Display, Error, Formatter}; + +/// Struct that represents a boundary. +/// This will be used to determine how the boundary conditions behaves when it interacts +/// with photon packets. +pub struct Boundary { + bounding_box: Cube, + top: BoundaryCondition, + bottom: BoundaryCondition, + north: BoundaryCondition, + east: BoundaryCondition, + south: BoundaryCondition, + west: BoundaryCondition, +} + +impl Boundary { + access!(bounding_box, bounding_box_mut: Cube); + access!(top, top_mut: BoundaryCondition); + access!(bottom, bottom_mut: BoundaryCondition); + access!(north, north_mut: BoundaryCondition); + access!(east, east_mut: BoundaryCondition); + access!(south, south_mut: BoundaryCondition); + access!(west, west_mut: BoundaryCondition); + + pub fn new(bounding_box: Cube) -> Self { + Self { + bounding_box, + top: BoundaryCondition::default(), + bottom: BoundaryCondition::default(), + north: BoundaryCondition::default(), + east: BoundaryCondition::default(), + south: BoundaryCondition::default(), + west: BoundaryCondition::default(), + } + } + + pub fn new_kill(bounding_box: Cube) -> Self { + Self { + bounding_box, + top: BoundaryCondition::Kill, + bottom: BoundaryCondition::Kill, + north: BoundaryCondition::Kill, + east: BoundaryCondition::Kill, + south: BoundaryCondition::Kill, + west: BoundaryCondition::Kill, + } + } + + pub fn new_reflect(bounding_box: Cube, reflect: Reflectance) -> Self { + Self { + bounding_box, + top: BoundaryCondition::Reflect(reflect.clone()), + bottom: BoundaryCondition::Reflect(reflect.clone()), + north: BoundaryCondition::Reflect(reflect.clone()), + east: BoundaryCondition::Reflect(reflect.clone()), + south: BoundaryCondition::Reflect(reflect.clone()), + west: BoundaryCondition::Reflect(reflect.clone()), + } + } + + pub fn new_periodic(bounding_box: Cube, padding: f64) -> Self { + Self { + bounding_box, + top: BoundaryCondition::Periodic(padding), + bottom: BoundaryCondition::Periodic(padding), + north: BoundaryCondition::Periodic(padding), + east: BoundaryCondition::Periodic(padding), + south: BoundaryCondition::Periodic(padding), + west: BoundaryCondition::Periodic(padding), + } + } + + #[inline] + pub fn apply<'a>(&self, rng: &mut ThreadRng, hit: &'a BoundaryHit<'a>, phot: &mut Photon) { + match hit.condition() { + BoundaryCondition::Kill => { + // Handle Kill variant + phot.kill(); + } + BoundaryCondition::Reflect(reflectance) => { + // Handle Reflect variant + + match reflectance.reflect(rng, &phot, &hit.get_hit()) { + Some(ray) => *phot.ray_mut() = ray, + None => phot.kill(), + } + } + #[cfg(not(feature = "mpi"))] + BoundaryCondition::Periodic(padding) => { + // Get the opposing boundary + self.set_ray_to_opposite_boundary(&mut phot.ray_mut(), hit.direction(), padding); + } + #[cfg(feature = "mpi")] + BoundaryCondition::Periodic(padding) => { + // Handle this variant in the case of MPI. + unimplemented!() + } + #[cfg(feature = "mpi")] + BoundaryCondition::MpiRank => { + // Handle MpiRank variant + } + }; + } + + /// Provides the translation to a Point3 (+/- a padding) to move it from one + /// boundary to the opposing boundary. The primary intended use of this code is + /// in the application of a periodic boundary on a single compute node. + #[inline] + pub fn get_periodic_translation(&self, bound: &BoundaryDirection, padding: &f64) -> Vec3 { + // First determine the vector component that we need to translate. + let trans_vec = match bound { + BoundaryDirection::Top | BoundaryDirection::Bottom => Vec3::new( + 0.0, + 0.0, + self.bounding_box.maxs()[2] - self.bounding_box.mins()[2] - padding, + ), + BoundaryDirection::North | BoundaryDirection::South => Vec3::new( + 0.0, + self.bounding_box.maxs()[1] - self.bounding_box.mins()[1] - padding, + 0.0, + ), + BoundaryDirection::East | BoundaryDirection::West => Vec3::new( + self.bounding_box.maxs()[0] - self.bounding_box.mins()[0] - padding, + 0.0, + 0.0, + ), + }; + + // Finally determine the direction of the resulting translation. + match bound { + BoundaryDirection::Top | BoundaryDirection::East | BoundaryDirection::North => { + -trans_vec + } + _ => trans_vec, + } + } + + #[inline] + pub fn set_ray_to_opposite_boundary( + &self, + ray: &mut Ray, + bound: &BoundaryDirection, + padding: &f64, + ) { + // First determine the axis that we need to translate. + let axis = match bound { + BoundaryDirection::Top | BoundaryDirection::Bottom => 2, + BoundaryDirection::North | BoundaryDirection::South => 1, + BoundaryDirection::East | BoundaryDirection::West => 0, + }; + + // Finally, set the position of the ray to the opposite boundary, +/- padding. + match bound { + BoundaryDirection::Top | BoundaryDirection::East | BoundaryDirection::North => { + ray.pos_mut()[axis] = self.bounding_box.mins()[axis] + padding; + } + _ => { + ray.pos_mut()[axis] = self.bounding_box.maxs()[axis] - padding; + } + } + } + + pub fn dist_boundary(&self, ray: &Ray) -> Option<(f64, BoundaryDirection)> { + if let Some((dist, side)) = self.bounding_box.dist_side(ray) { + debug_assert!(matches!(side, Side::Inside(_))); + + // Now we have to find the boundary that the ray is going to hit. + // We can do this by finding the max absolutel component value of the + // vector. Then, find the dir + return Some((dist, BoundaryDirection::ray_facing_boundary(ray))); + } + + None + } +} + +impl Display for Boundary { + fn fmt(&self, f: &mut Formatter<'_>) -> std::fmt::Result { + fmt_report!(f, self.bounding_box, "bounding box"); + fmt_report!(f, self.top, "top"); + fmt_report!(f, self.bottom, "bottom"); + fmt_report!(f, self.north, "north"); + fmt_report!(f, self.east, "east"); + fmt_report!(f, self.south, "south"); + fmt_report!(f, self.west, "west"); + Ok(()) + } +} + +/// Describing a boundary at which the action is triggered. +/// This will help determine how the boundary conditions behaves when it interacts +/// with photon packets. +#[derive(Debug, Eq, PartialEq, Clone)] +pub enum BoundaryDirection { + /// The boundary at the maximum z-value. + Top, + /// The boundary at the minimum z-value. + Bottom, + /// The boundary at the maximum y-value. + North, + /// The boundary at the maximum x-value. + East, + /// The boundary at the minimum y-value. + South, + /// The boundary at the minimum y-value. + West, +} + +impl BoundaryDirection { + pub fn opposing(&self) -> BoundaryDirection { + match self { + BoundaryDirection::Top => BoundaryDirection::Bottom, + BoundaryDirection::Bottom => BoundaryDirection::Top, + BoundaryDirection::North => BoundaryDirection::South, + BoundaryDirection::South => BoundaryDirection::North, + BoundaryDirection::East => BoundaryDirection::West, + BoundaryDirection::West => BoundaryDirection::East, + } + } + + pub fn normal_vector(&self) -> Dir3 { + match self { + BoundaryDirection::Top => Dir3::new(0.0, 0.0, 1.0), + BoundaryDirection::Bottom => Dir3::new(0.0, 0.0, -1.0), + BoundaryDirection::North => Dir3::new(0.0, 1.0, 0.0), + BoundaryDirection::South => Dir3::new(0.0, -1.0, 0.0), + BoundaryDirection::East => Dir3::new(1.0, 0.0, 0.0), + BoundaryDirection::West => Dir3::new(-1.0, 0.0, 0.0), + } + } + + // TODO: Do we need to consider what happens if we are oriented to an edge of corner? + + /// Determines the boundary which the ray is currently facing, and hence the + /// boundary which it is going to hit. + pub fn ray_facing_boundary(ray: &Ray) -> Self { + let direction = ray.dir(); + let abs_x = direction.x().abs(); + let abs_y = direction.y().abs(); + let abs_z = direction.z().abs(); + + if abs_x >= abs_y && abs_x >= abs_z { + if direction.x() > 0.0 { + BoundaryDirection::East + } else { + BoundaryDirection::West + } + } else if abs_y >= abs_x && abs_y >= abs_z { + if direction.y() > 0.0 { + BoundaryDirection::North + } else { + BoundaryDirection::South + } + } else { + if direction.z() > 0.0 { + BoundaryDirection::Top + } else { + BoundaryDirection::Bottom + } + } + } +} + +impl Display for BoundaryDirection { + fn fmt(&self, f: &mut Formatter) -> std::fmt::Result { + match self { + BoundaryDirection::Top => write!(f, "Top"), + BoundaryDirection::Bottom => write!(f, "Bottom"), + BoundaryDirection::North => write!(f, "North"), + BoundaryDirection::East => write!(f, "East"), + BoundaryDirection::South => write!(f, "South"), + BoundaryDirection::West => write!(f, "West"), + } + } +} + +#[derive(Default, Debug, PartialEq)] +pub enum BoundaryCondition { + /// Any photon packet that intersects with this boundary will be down-weighted + /// and removed from the simulation. + #[default] + Kill, + /// Any photon packet that intersects with this boundarty will be specularly + /// reflected back into the domain. + Reflect(Reflectance), + /// Any photon that intersects with this boundary will be transferred to the + /// opposing boundary and re-emitted + /// The number defines the padding distance from the oppising edge (to avoid instant re-collision). + Periodic(f64), + /// Photons that intersect this boundary will be collected, buffered and + /// transferred to the adjacent MPI rank. + #[cfg(feature = "mpi")] + MpiRank(usize), +} + + +impl Display for BoundaryCondition { + fn fmt(&self, fmt: &mut Formatter<'_>) -> std::fmt::Result { + match *self { + Self::Kill => { + writeln!(fmt, "Kill: ...")?; + Ok(()) + } + Self::Reflect(ref reflectance) => { + writeln!(fmt, "Reflector: ...")?; + fmt_report!(fmt, reflectance, "reflectance"); + Ok(()) + }, + Self::Periodic(padding) => { + writeln!(fmt, "Periodic: ...")?; + fmt_report!(fmt, padding, "padding"); + Ok(()) + }, + #[cfg(feature = "mpi")] + Self::MpiRank(rank) => { + writeln!(fmt, "MPI Rank Transfer: ...")?; + fmt_report!(fmt, padding, "destination rank"); + Ok(()) + } + } + } +} + +#[derive(Debug, PartialEq, Clone)] +pub struct BoundaryHit<'a> { + condition: &'a BoundaryCondition, + dist: f64, + direction: BoundaryDirection, +} + +impl<'a> BoundaryHit<'a> { + access!(condition: BoundaryCondition); + clone!(dist, dist_mut: f64); + access!(direction: BoundaryDirection); + + #[inline] + #[must_use] + pub fn new(condition: &'a BoundaryCondition, dist: f64, direction: BoundaryDirection) -> Self { + debug_assert!(dist > 0.0); + Self { + condition, + dist, + direction, + } + } + + pub fn get_hit(&self) -> Hit<'_, Attribute> { + Hit::new( + &Attribute::Mirror(0.0), + self.dist(), + Side::Inside(self.direction().normal_vector()), + ) + } +} + +impl<'a> Into>> for BoundaryHit<'a> { + fn into(self) -> Hit<'a, Attribute<'a>> { + // Not the most elegant implementation, as the tag is not used. + Hit::new( + &Attribute::Mirror(0.0), + self.dist(), + Side::Inside(self.direction().normal_vector()), + ) + } +} + +#[cfg(test)] +mod tests { + use super::*; + use crate::math::{Dir3, Point3}; + + #[test] + fn test_boundary_facing() { + let boundary = Boundary { + bounding_box: Cube::new(Point3::new(0.0, 0.0, 0.0), Point3::new(6.0, 8.0, 10.0)), + top: BoundaryCondition::default(), + bottom: BoundaryCondition::default(), + north: BoundaryCondition::default(), + east: BoundaryCondition::default(), + south: BoundaryCondition::default(), + west: BoundaryCondition::default(), + }; + + // Basic - Check a ray facing zenith. + let ray = Ray::new(Point3::new(5.0, 5.0, 5.0), Dir3::new(0.0, 0.0, 1.0)); + let (dist, bound) = boundary + .dist_boundary(&ray) + .expect("Ray not contained within domain. "); + + // Now do our first test. + assert_eq!(dist, 5.0); + assert_eq!(bound, BoundaryDirection::Top); + + // Update the ray to face nadir. + let ray = Ray::new(Point3::new(5.0, 5.0, 5.0), Dir3::new(0.0, 0.0, -1.0)); + let (dist, bound) = boundary + .dist_boundary(&ray) + .expect("Ray not contained within domain. "); + // Now test that we get the correct boundary. + assert_eq!(dist, 5.0); + assert_eq!(bound, BoundaryDirection::Bottom); + + // Now do a test for the east boundary. + // Update the ray to face east. + let ray = Ray::new(Point3::new(5.0, 5.0, 5.0), Dir3::new(1.0, 0.0, 0.0)); + let (dist, bound) = boundary + .dist_boundary(&ray) + .expect("Ray not contained within domain. "); + // Now test that we get the correct boundary. + assert_eq!(dist, 1.0); + assert_eq!(bound, BoundaryDirection::East); + + // Now do a test for the west boundary. + // Update the ray to face west. + let ray = Ray::new(Point3::new(5.0, 5.0, 5.0), Dir3::new(-1.0, 0.0, 0.0)); + let (dist, bound) = boundary + .dist_boundary(&ray) + .expect("Ray not contained within domain. "); + // Now test that we get the correct boundary. + assert_eq!(dist, 5.0); + assert_eq!(bound, BoundaryDirection::West); + + // Now do a test for the north boundary. + // Update the ray to face north. + let ray = Ray::new(Point3::new(5.0, 5.0, 5.0), Dir3::new(0.0, 1.0, 0.0)); + let (dist, bound) = boundary + .dist_boundary(&ray) + .expect("Ray not contained within domain. "); + // Now test that we get the correct boundary. + assert_eq!(dist, 3.0); + assert_eq!(bound, BoundaryDirection::North); + + // Now do a test for the south boundary. + // Update the ray to face south. + let ray = Ray::new(Point3::new(5.0, 5.0, 5.0), Dir3::new(0.0, -1.0, 0.0)); + let (dist, bound) = boundary + .dist_boundary(&ray) + .expect("Ray not contained within domain. "); + // Now test that we get the correct boundary. + assert_eq!(dist, 5.0); + assert_eq!(bound, BoundaryDirection::South); + + // Now let's do a test that checks an odd combination of directions. + let ray = Ray::new(Point3::new(5.0, 5.0, 5.0), Dir3::new(1.0, 3.0, 5.0)); + let (dist, bound) = boundary + .dist_boundary(&ray) + .expect("Ray not contained within domain. "); + // Now test that we get the correct boundary. + + assert_eq!( + dist, + (1.0_f64 + 3.0_f64 * 3.0_f64 + 5.0_f64 * 5.0_f64).sqrt() + ); + assert_eq!(bound, BoundaryDirection::Top); + } + + #[test] + fn test_periodic_boundary() { + let mut rng = rand::thread_rng(); + + // Setup a basic boundary to the simulation. + // Each side is a different length, and is periodic. + let boundary = Boundary { + bounding_box: Cube::new(Point3::new(0.0, 0.0, 0.0), Point3::new(6.0, 8.0, 10.0)), + top: BoundaryCondition::Periodic(0.0), + bottom: BoundaryCondition::Periodic(0.0), + north: BoundaryCondition::Periodic(0.0), + east: BoundaryCondition::Periodic(0.0), + south: BoundaryCondition::Periodic(0.0), + west: BoundaryCondition::Periodic(0.0), + }; + + // Test with padding of 0.0 + let incoming_ray = Ray::new(Point3::new(5.0, 5.0, 9.98), Dir3::new(0.0, 0.0, 1.0)); + let mut incoming_photon = Photon::new(incoming_ray, 550.0, 1.0); + + let (dist, bound) = boundary + .dist_boundary(incoming_photon.ray()) + .expect("Ray not contained within domain. "); + let bhit = BoundaryHit::new(&BoundaryCondition::Periodic(0.0), dist, bound); + boundary.apply(&mut rng, &bhit, &mut incoming_photon); + assert_eq!(*incoming_photon.ray().pos(), Point3::new(5.0, 5.0, 0.0)); + assert_eq!(*incoming_photon.ray().dir(), Dir3::new(0.0, 0.0, 1.0)); + + // Test with padding of 0.01 + let incoming_ray = Ray::new(Point3::new(5.0, 0.02, 5.0), Dir3::new(0.1, -0.9, 0.0)); + let mut incoming_photon = Photon::new(incoming_ray, 550.0, 1.0); + + let (dist, bound) = boundary + .dist_boundary(incoming_photon.ray()) + .expect("Ray not contained within domain. "); + let bhit = BoundaryHit::new(&BoundaryCondition::Periodic(0.01), dist, bound); + boundary.apply(&mut rng, &bhit, &mut incoming_photon); + assert_eq!(*incoming_photon.ray().pos(), Point3::new(5.0, 7.99, 5.0)); + assert_eq!(*incoming_photon.ray().dir(), Dir3::new(0.1, -0.9, 0.0)); + } +} diff --git a/src/geom/domain/mod.rs b/src/geom/domain/mod.rs index 1a65c64..3d74bb7 100644 --- a/src/geom/domain/mod.rs +++ b/src/geom/domain/mod.rs @@ -1,5 +1,6 @@ //! Domain module. +pub mod boundary; pub mod grid; pub mod grid_builder; pub mod surface; @@ -9,6 +10,6 @@ pub mod tree; pub mod tree_settings; pub use self::{ - grid::*, grid_builder::*, surface::*, surface_linker::*, surface_linker_loader::*, tree::*, - tree_settings::*, + boundary::*, grid::*, grid_builder::*, surface::*, surface_linker::*, + surface_linker_loader::*, tree::*, tree_settings::*, }; diff --git a/src/sim/engine/engines/fluorescence.rs b/src/sim/engine/engines/fluorescence.rs index 0b00431..5c16abe 100644 --- a/src/sim/engine/engines/fluorescence.rs +++ b/src/sim/engine/engines/fluorescence.rs @@ -76,9 +76,10 @@ pub fn fluorescence( let surf_hit = input .tree .scan(phot.ray().clone(), bump_dist, voxel_dist.min(scat_dist)); + let boundary_hit = None; // Event handling. - match Event::new(voxel_dist, scat_dist, surf_hit, bump_dist) { + match Event::new(voxel_dist, scat_dist, surf_hit, boundary_hit, bump_dist) { Event::Voxel(dist) => travel(&mut data, &mut phot, &env, index, dist + bump_dist), Event::Scattering(dist) => { travel(&mut data, &mut phot, &env, index, dist); @@ -88,6 +89,10 @@ pub fn fluorescence( travel(&mut data, &mut phot, &env, index, hit.dist()); surface(&mut rng, &hit, &mut phot, &mut local, &mut data); travel(&mut data, &mut phot, &env, index, bump_dist); + }, + Event::Boundary(boundary_hit) => { + travel(&mut data, &mut phot, &env, index, boundary_hit.dist()); + input.bound.apply(rng, &boundary_hit, &mut phot); } } diff --git a/src/sim/engine/engines/photo.rs b/src/sim/engine/engines/photo.rs index 8fb9d34..982118b 100644 --- a/src/sim/engine/engines/photo.rs +++ b/src/sim/engine/engines/photo.rs @@ -67,9 +67,10 @@ pub fn photo( let surf_hit = input .tree .scan(phot.ray().clone(), bump_dist, voxel_dist.min(scat_dist)); + let boundary_hit = None; // Event handling. - match Event::new(voxel_dist, scat_dist, surf_hit, bump_dist) { + match Event::new(voxel_dist, scat_dist, surf_hit, boundary_hit, bump_dist) { Event::Voxel(dist) => travel(&mut data, &mut phot, &env, index, dist + bump_dist), Event::Scattering(dist) => { travel(&mut data, &mut phot, &env, index, dist); @@ -95,6 +96,10 @@ pub fn photo( travel(&mut data, &mut phot, &env, index, hit.dist()); surface(&mut rng, &hit, &mut phot, &mut env, &mut data); travel(&mut data, &mut phot, &env, index, bump_dist); + }, + Event::Boundary(boundary_hit) => { + travel(&mut data, &mut phot, &env, index, boundary_hit.dist()); + input.bound.apply(rng, &boundary_hit, &mut phot); } } diff --git a/src/sim/engine/engines/raman.rs b/src/sim/engine/engines/raman.rs index 739a065..83bd6a0 100644 --- a/src/sim/engine/engines/raman.rs +++ b/src/sim/engine/engines/raman.rs @@ -64,9 +64,10 @@ pub fn raman( let surf_hit = input .tree .scan(phot.ray().clone(), bump_dist, voxel_dist.min(scat_dist)); + let boundary_hit = None; // Event handling. - match Event::new(voxel_dist, scat_dist, surf_hit, bump_dist) { + match Event::new(voxel_dist, scat_dist, surf_hit, boundary_hit, bump_dist) { Event::Voxel(dist) => travel(&mut data, &mut phot, &env, index, dist + bump_dist), Event::Scattering(dist) => { travel(&mut data, &mut phot, &env, index, dist); @@ -84,6 +85,10 @@ pub fn raman( travel(&mut data, &mut phot, &env, index, hit.dist()); surface(&mut rng, &hit, &mut phot, &mut env, &mut data); travel(&mut data, &mut phot, &env, index, bump_dist); + }, + Event::Boundary(boundary_hit) => { + travel(&mut data, &mut phot, &env, index, boundary_hit.dist()); + input.bound.apply(rng, &boundary_hit, &mut phot); } } diff --git a/src/sim/engine/engines/standard.rs b/src/sim/engine/engines/standard.rs index 63a837b..bc93620 100644 --- a/src/sim/engine/engines/standard.rs +++ b/src/sim/engine/engines/standard.rs @@ -1,7 +1,7 @@ //! Standard photon-lifetime engine function. use crate::{ - geom::Trace, + geom::{Boundary, Trace}, phys::Photon, sim::{scatter::scatter, surface::surface, travel::travel, Event, Input, Output}, }; @@ -56,9 +56,10 @@ pub fn standard(input: &Input, mut data: &mut Output, mut rng: &mut ThreadRng, m let surf_hit = input .tree .scan(phot.ray().clone(), bump_dist, voxel_dist.min(scat_dist)); + let boundary_hit = None; // Event handling. - match Event::new(voxel_dist, scat_dist, surf_hit, bump_dist) { + match Event::new(voxel_dist, scat_dist, surf_hit, boundary_hit, bump_dist) { Event::Voxel(dist) => travel(&mut data, &mut phot, &env, index, dist + bump_dist), Event::Scattering(dist) => { travel(&mut data, &mut phot, &env, index, dist); @@ -68,6 +69,10 @@ pub fn standard(input: &Input, mut data: &mut Output, mut rng: &mut ThreadRng, m travel(&mut data, &mut phot, &env, index, hit.dist()); surface(&mut rng, &hit, &mut phot, &mut env, &mut data); travel(&mut data, &mut phot, &env, index, bump_dist); + }, + Event::Boundary(boundary_hit) => { + travel(&mut data, &mut phot, &env, index, boundary_hit.dist()); + input.bound.apply(rng, &boundary_hit, &mut phot); } } diff --git a/src/sim/event.rs b/src/sim/event.rs index b1f28b4..619d4d4 100644 --- a/src/sim/event.rs +++ b/src/sim/event.rs @@ -1,6 +1,6 @@ //! Event enumeration. -use crate::geom::Hit; +use crate::geom::{boundary, BoundaryHit, Hit}; /// Event determination enumeration. #[derive(PartialEq, Debug)] @@ -12,6 +12,8 @@ pub enum Event<'a, T> Scattering(f64), /// Surface hit. Surface(Hit<'a, T>), + /// Boundary + Boundary(BoundaryHit<'a>) } impl<'a, T> Event<'a, T> { @@ -23,6 +25,7 @@ impl<'a, T> Event<'a, T> { voxel_dist: f64, scat_dist: f64, surf_hit: Option>, + boundary_hit: Option>, bump_dist: f64, ) -> Self { debug_assert!(voxel_dist > 0.0); @@ -43,6 +46,13 @@ impl<'a, T> Event<'a, T> { return Self::Surface(hit); } + // We should be able to safely assume that if there were geometry in the + if let Some(bhit) = boundary_hit { + if bhit.dist() < scat_dist && bhit.dist() < (voxel_dist + bump_dist) { + return Self::Boundary(bhit) + } + } + if scat_dist < (voxel_dist + bump_dist) { return Self::Scattering(scat_dist); } @@ -54,7 +64,7 @@ impl<'a, T> Event<'a, T> { mod tests { use crate::{ sim::Attribute, - geom::Side, + geom::{Side, BoundaryCondition}, math::Dir3, }; use super::*; @@ -63,7 +73,7 @@ mod tests { #[test] fn test_new_surface_hit() { let surf_hit = Some(Hit::new(&Attribute::Mirror(0.5), 1.0, Side::Outside(Dir3::new(1.0, 0.0, 0.0)))); - let event = Event::new(2.0, 3.0, surf_hit, 0.5); + let event = Event::new(2.0, 3.0, surf_hit, None, 0.5); // Check each of the components of the event. if let Event::Surface(hit) = event { @@ -77,14 +87,22 @@ mod tests { #[test] fn test_new_voxel_collision() { - let event: Event<'_, Attribute> = Event::new(2.0, 3.0, None, 0.5); + let event: Event<'_, Attribute> = Event::new(2.0, 3.0, None, None, 0.5); assert_eq!(event, Event::Voxel(2.0)); } #[test] fn test_new_scattering_event() { let surf_hit = Some(Hit::new(&Attribute::Mirror(0.5), 2.0, Side::Outside(Dir3::new(1.0, 0.0, 0.0)))); - let event = Event::new(2.0, 1.0, surf_hit, 4.0); + let event = Event::new(2.0, 1.0, surf_hit, None, 0.5); assert_eq!(event, Event::Scattering(1.0)); } + + #[test] + fn test_new_boundary_event() { + + let bhit = BoundaryHit::new(&BoundaryCondition::Periodic(0.0), 0.1, boundary::BoundaryDirection::North); + let event: Event<'_, Attribute<'_>> = Event::new(2.0, 1.0, None, Some(bhit.clone()), 0.5); + assert_eq!(event, Event::Boundary(bhit)); + } } diff --git a/src/sim/input.rs b/src/sim/input.rs index 4b4b02a..0774d0f 100644 --- a/src/sim/input.rs +++ b/src/sim/input.rs @@ -2,7 +2,7 @@ use crate::{ fmt_report, - geom::{Grid, Tree}, + geom::{Boundary, Grid, Tree}, ord::{Register, Set}, phys::{Light, Material}, sim::{Attribute, Settings}, @@ -26,6 +26,8 @@ pub struct Input<'a> { pub grid: &'a Grid, /// General settings. pub sett: &'a Settings, + /// Boundary for the simulation. + pub bound: &'a Boundary, } impl<'a> Input<'a> { @@ -40,6 +42,7 @@ impl<'a> Input<'a> { tree: &'a Tree, grid: &'a Grid, sett: &'a Settings, + bound: &'a Boundary ) -> Self { Self { spec_reg, @@ -49,6 +52,7 @@ impl<'a> Input<'a> { tree, grid, sett, + bound, } } } @@ -64,6 +68,7 @@ impl Display for Input<'_> { fmt_report!(fmt, self.tree, "hit-scan tree"); fmt_report!(fmt, self.grid, "measurement grid"); fmt_report!(fmt, self.sett, "settings"); + fmt_report!(fmt, self.bound, "boundary"); Ok(()) } } From 1b7921dc0c25513534fa0630671ddc38e50a76b7 Mon Sep 17 00:00:00 2001 From: Sam Morrell Date: Tue, 16 Apr 2024 01:05:54 +0100 Subject: [PATCH 02/16] Initial integration of boundary condition into MCRT - Improvement of boundary condition structs to work with MCRT. - Fixed bug where boundary hit finding would not be tolerant of rounding. - Integration of boundary condition into engines and MCRT main. - Added supporting ray-plane intersect function to support changes. - Removed dependence of engine and travel upon output grid that stretches across whole domain. Making way for more flexible output in the future. --- src/bin/mcrt.rs | 3 +- src/geom/domain/boundary.rs | 225 +++++++++++++++++-------- src/geom/mod.rs | 1 + src/geom/plane.rs | 21 +++ src/sim/engine/engines/fluorescence.rs | 26 +-- src/sim/engine/engines/photo.rs | 24 +-- src/sim/engine/engines/raman.rs | 24 +-- src/sim/engine/engines/standard.rs | 30 ++-- src/sim/event.rs | 31 ++-- src/sim/mod.rs | 3 +- src/sim/travel.rs | 16 +- 11 files changed, 267 insertions(+), 137 deletions(-) create mode 100644 src/geom/plane.rs diff --git a/src/bin/mcrt.rs b/src/bin/mcrt.rs index 8cb7394..72d92bb 100644 --- a/src/bin/mcrt.rs +++ b/src/bin/mcrt.rs @@ -10,9 +10,8 @@ use aetherus::{ args, data::Histogram, fs::{File, Load, Save}, - geom::{Grid, Tree, Boundary, BoundaryCondition}, + geom::{Grid, Tree, Boundary}, img::{Colour, Image}, - math::{Point3}, ord::{Build, Link, Register, Set, cartesian::{X, Y}}, report, sim::{ diff --git a/src/geom/domain/boundary.rs b/src/geom/domain/boundary.rs index 328aba3..8eaca6c 100644 --- a/src/geom/domain/boundary.rs +++ b/src/geom/domain/boundary.rs @@ -1,8 +1,8 @@ use crate::{ - access, clone, fmt_report, geom::{Cube, Hit, Ray, Side, Trace}, math::{Dir3, Vec3}, phys::{Photon, Reflectance}, sim::Attribute + access, clone, fmt_report, geom::{plane::ray_plane_intersection, Cube, Hit, Ray, Side, Trace}, math::{Dir3, Point3, Vec3}, phys::{Photon, Reflectance}, sim::Attribute }; use rand::rngs::ThreadRng; -use std::fmt::{Display, Error, Formatter}; +use std::fmt::{Display, Formatter}; /// Struct that represents a boundary. /// This will be used to determine how the boundary conditions behaves when it interacts @@ -83,7 +83,6 @@ impl Boundary { } BoundaryCondition::Reflect(reflectance) => { // Handle Reflect variant - match reflectance.reflect(rng, &phot, &hit.get_hit()) { Some(ray) => *phot.ray_mut() = ray, None => phot.kill(), @@ -164,18 +163,67 @@ impl Boundary { } } - pub fn dist_boundary(&self, ray: &Ray) -> Option<(f64, BoundaryDirection)> { - if let Some((dist, side)) = self.bounding_box.dist_side(ray) { - debug_assert!(matches!(side, Side::Inside(_))); - + pub fn dist_boundary(&self, ray: &Ray) -> Option { + if let Some((dist, _side)) = self.bounding_box.dist_side(ray) { // Now we have to find the boundary that the ray is going to hit. // We can do this by finding the max absolutel component value of the // vector. Then, find the dir - return Some((dist, BoundaryDirection::ray_facing_boundary(ray))); + let dir = self.boundary_direction(&ray).expect("Ray outside of boundary. "); + return Some(BoundaryHit::new( + self.condition_for_boundary(&dir), + dist, + dir, + )) } None } + + #[inline] + pub fn condition_for_boundary(&self, dir: &BoundaryDirection) -> &BoundaryCondition { + match &dir { + BoundaryDirection::Bottom => self.bottom(), + BoundaryDirection::Top => self.top(), + BoundaryDirection::North => self.north(), + BoundaryDirection::South => self.south(), + BoundaryDirection::West => self.west(), + BoundaryDirection::East => self.east(), + } + } + + #[inline] + pub fn contains(&self, p: &Point3) -> bool { + self.bounding_box.contains(p) + } + + pub fn boundary_direction(&self, ray: &Ray) -> Option { + const TOLLERANCE: f64 = 1E-10; + let dirs = [ + BoundaryDirection::Bottom, + BoundaryDirection::Top, + BoundaryDirection::South, + BoundaryDirection::North, + BoundaryDirection::West, + BoundaryDirection::East, + ]; + + let min_point = Point3::new(self.bounding_box().mins()[0], self.bounding_box().mins()[1], self.bounding_box().mins()[2]); + let max_point = Point3::new(self.bounding_box().maxs()[0], self.bounding_box().maxs()[1], self.bounding_box().maxs()[2]); + let mut facing_dir = None; + for (i, &dir) in dirs.iter().enumerate() { + if let Some(point) = ray_plane_intersection(&ray, if i % 2 == 0 { min_point } else { max_point }, dir.normal_vector()) { + if point.x() >= min_point.x() - TOLLERANCE && point.x() <= max_point.x() + TOLLERANCE && + point.y() >= min_point.y() - TOLLERANCE && point.y() <= max_point.y() + TOLLERANCE && + point.z() >= min_point.z() - TOLLERANCE && point.z() <= max_point.z() + TOLLERANCE + { + facing_dir = Some(dir.clone()); + break; + } + } + } + + return facing_dir; + } } impl Display for Boundary { @@ -194,7 +242,7 @@ impl Display for Boundary { /// Describing a boundary at which the action is triggered. /// This will help determine how the boundary conditions behaves when it interacts /// with photon packets. -#[derive(Debug, Eq, PartialEq, Clone)] +#[derive(Debug, Eq, PartialEq, Clone, Copy)] pub enum BoundaryDirection { /// The boundary at the maximum z-value. Top, @@ -224,43 +272,12 @@ impl BoundaryDirection { pub fn normal_vector(&self) -> Dir3 { match self { - BoundaryDirection::Top => Dir3::new(0.0, 0.0, 1.0), - BoundaryDirection::Bottom => Dir3::new(0.0, 0.0, -1.0), - BoundaryDirection::North => Dir3::new(0.0, 1.0, 0.0), - BoundaryDirection::South => Dir3::new(0.0, -1.0, 0.0), - BoundaryDirection::East => Dir3::new(1.0, 0.0, 0.0), - BoundaryDirection::West => Dir3::new(-1.0, 0.0, 0.0), - } - } - - // TODO: Do we need to consider what happens if we are oriented to an edge of corner? - - /// Determines the boundary which the ray is currently facing, and hence the - /// boundary which it is going to hit. - pub fn ray_facing_boundary(ray: &Ray) -> Self { - let direction = ray.dir(); - let abs_x = direction.x().abs(); - let abs_y = direction.y().abs(); - let abs_z = direction.z().abs(); - - if abs_x >= abs_y && abs_x >= abs_z { - if direction.x() > 0.0 { - BoundaryDirection::East - } else { - BoundaryDirection::West - } - } else if abs_y >= abs_x && abs_y >= abs_z { - if direction.y() > 0.0 { - BoundaryDirection::North - } else { - BoundaryDirection::South - } - } else { - if direction.z() > 0.0 { - BoundaryDirection::Top - } else { - BoundaryDirection::Bottom - } + BoundaryDirection::Top => Dir3::new(0.0, 0.0, -1.0), + BoundaryDirection::Bottom => Dir3::new(0.0, 0.0, 1.0), + BoundaryDirection::North => Dir3::new(0.0, -1.0, 0.0), + BoundaryDirection::South => Dir3::new(0.0, 1.0, 0.0), + BoundaryDirection::East => Dir3::new(-1.0, 0.0, 0.0), + BoundaryDirection::West => Dir3::new(1.0, 0.0, 0.0), } } } @@ -387,75 +404,75 @@ mod tests { // Basic - Check a ray facing zenith. let ray = Ray::new(Point3::new(5.0, 5.0, 5.0), Dir3::new(0.0, 0.0, 1.0)); - let (dist, bound) = boundary + let boundary_hit = boundary .dist_boundary(&ray) .expect("Ray not contained within domain. "); // Now do our first test. - assert_eq!(dist, 5.0); - assert_eq!(bound, BoundaryDirection::Top); + assert_eq!(boundary_hit.dist(), 5.0); + assert_eq!(*boundary_hit.direction(), BoundaryDirection::Top); // Update the ray to face nadir. let ray = Ray::new(Point3::new(5.0, 5.0, 5.0), Dir3::new(0.0, 0.0, -1.0)); - let (dist, bound) = boundary + let boundary_hit = boundary .dist_boundary(&ray) .expect("Ray not contained within domain. "); // Now test that we get the correct boundary. - assert_eq!(dist, 5.0); - assert_eq!(bound, BoundaryDirection::Bottom); + assert_eq!(boundary_hit.dist(), 5.0); + assert_eq!(*boundary_hit.direction(), BoundaryDirection::Bottom); // Now do a test for the east boundary. // Update the ray to face east. let ray = Ray::new(Point3::new(5.0, 5.0, 5.0), Dir3::new(1.0, 0.0, 0.0)); - let (dist, bound) = boundary + let boundary_hit = boundary .dist_boundary(&ray) .expect("Ray not contained within domain. "); // Now test that we get the correct boundary. - assert_eq!(dist, 1.0); - assert_eq!(bound, BoundaryDirection::East); + assert_eq!(boundary_hit.dist(), 1.0); + assert_eq!(*boundary_hit.direction(), BoundaryDirection::East); // Now do a test for the west boundary. // Update the ray to face west. let ray = Ray::new(Point3::new(5.0, 5.0, 5.0), Dir3::new(-1.0, 0.0, 0.0)); - let (dist, bound) = boundary + let boundary_hit = boundary .dist_boundary(&ray) .expect("Ray not contained within domain. "); // Now test that we get the correct boundary. - assert_eq!(dist, 5.0); - assert_eq!(bound, BoundaryDirection::West); + assert_eq!(boundary_hit.dist(), 5.0); + assert_eq!(*boundary_hit.direction(), BoundaryDirection::West); // Now do a test for the north boundary. // Update the ray to face north. let ray = Ray::new(Point3::new(5.0, 5.0, 5.0), Dir3::new(0.0, 1.0, 0.0)); - let (dist, bound) = boundary + let boundary_hit = boundary .dist_boundary(&ray) .expect("Ray not contained within domain. "); // Now test that we get the correct boundary. - assert_eq!(dist, 3.0); - assert_eq!(bound, BoundaryDirection::North); + assert_eq!(boundary_hit.dist(), 3.0); + assert_eq!(*boundary_hit.direction(), BoundaryDirection::North); // Now do a test for the south boundary. // Update the ray to face south. let ray = Ray::new(Point3::new(5.0, 5.0, 5.0), Dir3::new(0.0, -1.0, 0.0)); - let (dist, bound) = boundary + let boundary_hit = boundary .dist_boundary(&ray) .expect("Ray not contained within domain. "); // Now test that we get the correct boundary. - assert_eq!(dist, 5.0); - assert_eq!(bound, BoundaryDirection::South); + assert_eq!(boundary_hit.dist(), 5.0); + assert_eq!(*boundary_hit.direction(), BoundaryDirection::South); // Now let's do a test that checks an odd combination of directions. let ray = Ray::new(Point3::new(5.0, 5.0, 5.0), Dir3::new(1.0, 3.0, 5.0)); - let (dist, bound) = boundary + let boundary_hit = boundary .dist_boundary(&ray) .expect("Ray not contained within domain. "); // Now test that we get the correct boundary. assert_eq!( - dist, + boundary_hit.dist(), (1.0_f64 + 3.0_f64 * 3.0_f64 + 5.0_f64 * 5.0_f64).sqrt() ); - assert_eq!(bound, BoundaryDirection::Top); + assert_eq!(*boundary_hit.direction(), BoundaryDirection::Top); } #[test] @@ -478,24 +495,90 @@ mod tests { let incoming_ray = Ray::new(Point3::new(5.0, 5.0, 9.98), Dir3::new(0.0, 0.0, 1.0)); let mut incoming_photon = Photon::new(incoming_ray, 550.0, 1.0); - let (dist, bound) = boundary + let bhit = boundary .dist_boundary(incoming_photon.ray()) .expect("Ray not contained within domain. "); - let bhit = BoundaryHit::new(&BoundaryCondition::Periodic(0.0), dist, bound); boundary.apply(&mut rng, &bhit, &mut incoming_photon); assert_eq!(*incoming_photon.ray().pos(), Point3::new(5.0, 5.0, 0.0)); assert_eq!(*incoming_photon.ray().dir(), Dir3::new(0.0, 0.0, 1.0)); + let boundary = Boundary { + bounding_box: Cube::new(Point3::new(0.0, 0.0, 0.0), Point3::new(6.0, 8.0, 10.0)), + top: BoundaryCondition::Periodic(0.01), + bottom: BoundaryCondition::Periodic(0.01), + north: BoundaryCondition::Periodic(0.01), + east: BoundaryCondition::Periodic(0.01), + south: BoundaryCondition::Periodic(0.01), + west: BoundaryCondition::Periodic(0.01), + }; + // Test with padding of 0.01 let incoming_ray = Ray::new(Point3::new(5.0, 0.02, 5.0), Dir3::new(0.1, -0.9, 0.0)); let mut incoming_photon = Photon::new(incoming_ray, 550.0, 1.0); - let (dist, bound) = boundary + let bhit = boundary .dist_boundary(incoming_photon.ray()) .expect("Ray not contained within domain. "); - let bhit = BoundaryHit::new(&BoundaryCondition::Periodic(0.01), dist, bound); boundary.apply(&mut rng, &bhit, &mut incoming_photon); assert_eq!(*incoming_photon.ray().pos(), Point3::new(5.0, 7.99, 5.0)); assert_eq!(*incoming_photon.ray().dir(), Dir3::new(0.1, -0.9, 0.0)); } + + #[test] + fn test_cube_boundary_direction() { + // Setup a basic boundary to the simulation. + // Each side is a different length, and is periodic. + let boundary = Boundary { + bounding_box: Cube::new(Point3::new(0.0, 0.0, 0.0), Point3::new(1.0, 1.0, 1.0)), + top: BoundaryCondition::Periodic(0.0), + bottom: BoundaryCondition::Periodic(0.0), + north: BoundaryCondition::Periodic(0.0), + east: BoundaryCondition::Periodic(0.0), + south: BoundaryCondition::Periodic(0.0), + west: BoundaryCondition::Periodic(0.0), + }; + + let incoming_ray = Ray::new(Point3::new(0.5, 0.5, 0.5), Dir3::new(1.0, 0.0, 0.0)); + let b = boundary.boundary_direction(&incoming_ray); + assert_eq!(b, Some(BoundaryDirection::East)); + + let incoming_ray = Ray::new(Point3::new(0.5, 0.5, 0.5), Dir3::new(-1.0, 0.0, 0.0)); + let b = boundary.boundary_direction(&incoming_ray); + assert_eq!(b, Some(BoundaryDirection::West)); + + let incoming_ray = Ray::new(Point3::new(0.5, 0.5, 0.5), Dir3::new(0.0, 1.0, 0.0)); + let b = boundary.boundary_direction(&incoming_ray); + assert_eq!(b, Some(BoundaryDirection::North)); + + let incoming_ray = Ray::new(Point3::new(0.5, 0.5, 0.5), Dir3::new(0.0, -1.0, 0.0)); + let b = boundary.boundary_direction(&incoming_ray); + assert_eq!(b, Some(BoundaryDirection::South)); + + let incoming_ray = Ray::new(Point3::new(0.5, 0.5, 0.5), Dir3::new(0.0, 0.0, 1.0)); + let b = boundary.boundary_direction(&incoming_ray); + assert_eq!(b, Some(BoundaryDirection::Top)); + + let incoming_ray = Ray::new(Point3::new(0.5, 0.5, 0.5), Dir3::new(0.0, 0.0, -1.0)); + let b = boundary.boundary_direction(&incoming_ray); + assert_eq!(b, Some(BoundaryDirection::Bottom)); + } + + #[test] + fn test_boundary_edge() { + // Setup a basic boundary to the simulation. + // Each side is a different length, and is periodic. + let boundary = Boundary { + bounding_box: Cube::new(Point3::new(0.0, 0.0, 0.0), Point3::new(1.0, 1.0, 1.0)), + top: BoundaryCondition::Periodic(0.0), + bottom: BoundaryCondition::Periodic(0.0), + north: BoundaryCondition::Periodic(0.0), + east: BoundaryCondition::Periodic(0.0), + south: BoundaryCondition::Periodic(0.0), + west: BoundaryCondition::Periodic(0.0), + }; + + let incoming_ray = Ray::new(Point3::new(1.0, 0.5, 0.5), Dir3::new(1.0, 0.0, 0.0)); + let b = boundary.boundary_direction(&incoming_ray); + assert_eq!(b, Some(BoundaryDirection::East)); + } } diff --git a/src/geom/mod.rs b/src/geom/mod.rs index 5a86e95..d3a53f7 100644 --- a/src/geom/mod.rs +++ b/src/geom/mod.rs @@ -1,5 +1,6 @@ pub mod cast; pub mod domain; +pub mod plane; pub mod properties; pub mod rt; pub mod shape; diff --git a/src/geom/plane.rs b/src/geom/plane.rs new file mode 100644 index 0000000..b24013f --- /dev/null +++ b/src/geom/plane.rs @@ -0,0 +1,21 @@ +use crate::{ + geom::Ray, + math::{Dir3, Point3} +}; + +pub fn ray_plane_intersection(ray: &Ray, plane_point: Point3, plane_normal: Dir3) -> Option { + let p0_q0 = ray.pos() - plane_point; + let denom = ray.dir().dot(&plane_normal); + if denom.abs() < 1e-6 { + // The ray and the plane are parallel -- no intersection. + return None; + } + let dist = -p0_q0.dot(&plane_normal.into()) / denom; + if dist < 0.0 { + // Intersection is behind the ray, in the opposite direction to the dir vector. + return None; + } + + let ray_dir_vec = Point3::new(ray.dir().x(), ray.dir().y(), ray.dir().z()); + Some(ray_dir_vec * dist + *ray.pos()) +} \ No newline at end of file diff --git a/src/sim/engine/engines/fluorescence.rs b/src/sim/engine/engines/fluorescence.rs index 5c16abe..b571193 100644 --- a/src/sim/engine/engines/fluorescence.rs +++ b/src/sim/engine/engines/fluorescence.rs @@ -42,7 +42,7 @@ pub fn fluorescence( // Main event loop. let mut num_loops = 0; - while let Some((index, voxel)) = input.grid.gen_index_voxel(phot.ray().pos()) { + while input.bound.contains(phot.ray().pos()) { // Loop limit check. if num_loops >= loop_limit { println!("[WARN] : Terminating photon: loop limit reached."); @@ -60,38 +60,42 @@ pub fn fluorescence( } // Local variable modifications. + let index = input.grid.gen_index_voxel(phot.ray().pos()); env = Local::new( local.ref_index(), local.scat_coeff(), local.abs_coeff(), - mu_shift.mul_add(flu_concs[index], local.shift_coeff()), + mu_shift.mul_add(flu_concs[index.as_ref().unwrap().0], local.shift_coeff()), local.asym(), ); // Interaction distances. - let voxel_dist = voxel - .dist(phot.ray()) - .expect("Could not determine voxel distance."); + let voxel_dist = match &index { + Some((_index, voxel)) => { + voxel.dist(phot.ray()).expect("Could not determine voxel distance.") + }, + None => f64::INFINITY, + }; let scat_dist = -(rng.gen::()).ln() / env.inter_coeff(); let surf_hit = input .tree .scan(phot.ray().clone(), bump_dist, voxel_dist.min(scat_dist)); - let boundary_hit = None; + let boundary_hit = input.bound.dist_boundary(phot.ray()).expect("Photon not contained in boundary. "); // Event handling. match Event::new(voxel_dist, scat_dist, surf_hit, boundary_hit, bump_dist) { - Event::Voxel(dist) => travel(&mut data, &mut phot, &env, index, dist + bump_dist), + Event::Voxel(dist) => travel(&mut data, &mut phot, &env, &index, dist + bump_dist), Event::Scattering(dist) => { - travel(&mut data, &mut phot, &env, index, dist); + travel(&mut data, &mut phot, &env, &index, dist); scatter(&mut rng, &mut phot, &env); } Event::Surface(hit) => { - travel(&mut data, &mut phot, &env, index, hit.dist()); + travel(&mut data, &mut phot, &env, &index, hit.dist()); surface(&mut rng, &hit, &mut phot, &mut local, &mut data); - travel(&mut data, &mut phot, &env, index, bump_dist); + travel(&mut data, &mut phot, &env, &index, bump_dist); }, Event::Boundary(boundary_hit) => { - travel(&mut data, &mut phot, &env, index, boundary_hit.dist()); + travel(&mut data, &mut phot, &env, &index, boundary_hit.dist()); input.bound.apply(rng, &boundary_hit, &mut phot); } } diff --git a/src/sim/engine/engines/photo.rs b/src/sim/engine/engines/photo.rs index 982118b..2cc3014 100644 --- a/src/sim/engine/engines/photo.rs +++ b/src/sim/engine/engines/photo.rs @@ -42,7 +42,7 @@ pub fn photo( // Main event loop. let mut num_loops = 0; - while let Some((index, voxel)) = input.grid.gen_index_voxel(phot.ray().pos()) { + while input.bound.contains(phot.ray().pos()) { // Loop limit check. if num_loops >= loop_limit { println!("[WARN] : Terminating photon: loop limit reached."); @@ -60,20 +60,24 @@ pub fn photo( } // Interaction distances. - let voxel_dist = voxel - .dist(phot.ray()) - .expect("Could not determine voxel distance."); + let index = input.grid.gen_index_voxel(phot.ray().pos()); + let voxel_dist = match &index { + Some((_index, voxel)) => { + voxel.dist(phot.ray()).expect("Could not determine voxel distance.") + }, + None => f64::INFINITY, + }; let scat_dist = -(rng.gen::()).ln() / env.inter_coeff(); let surf_hit = input .tree .scan(phot.ray().clone(), bump_dist, voxel_dist.min(scat_dist)); - let boundary_hit = None; + let boundary_hit = input.bound.dist_boundary(phot.ray()).expect("Photon not contained in boundary. "); // Event handling. match Event::new(voxel_dist, scat_dist, surf_hit, boundary_hit, bump_dist) { - Event::Voxel(dist) => travel(&mut data, &mut phot, &env, index, dist + bump_dist), + Event::Voxel(dist) => travel(&mut data, &mut phot, &env, &index, dist + bump_dist), Event::Scattering(dist) => { - travel(&mut data, &mut phot, &env, index, dist); + travel(&mut data, &mut phot, &env, &index, dist); // Capture. for (frame, photo) in frames.iter().zip(data.photos.iter_mut()) { @@ -93,12 +97,12 @@ pub fn photo( scatter(&mut rng, &mut phot, &env); } Event::Surface(hit) => { - travel(&mut data, &mut phot, &env, index, hit.dist()); + travel(&mut data, &mut phot, &env, &index, hit.dist()); surface(&mut rng, &hit, &mut phot, &mut env, &mut data); - travel(&mut data, &mut phot, &env, index, bump_dist); + travel(&mut data, &mut phot, &env, &index, bump_dist); }, Event::Boundary(boundary_hit) => { - travel(&mut data, &mut phot, &env, index, boundary_hit.dist()); + travel(&mut data, &mut phot, &env, &index, boundary_hit.dist()); input.bound.apply(rng, &boundary_hit, &mut phot); } } diff --git a/src/sim/engine/engines/raman.rs b/src/sim/engine/engines/raman.rs index 83bd6a0..de1431f 100644 --- a/src/sim/engine/engines/raman.rs +++ b/src/sim/engine/engines/raman.rs @@ -39,7 +39,7 @@ pub fn raman( // Main event loop. let mut num_loops = 0; - while let Some((index, voxel)) = input.grid.gen_index_voxel(phot.ray().pos()) { + while input.bound.contains(phot.ray().pos()) { // Loop limit check. if num_loops >= loop_limit { println!("[WARN] : Terminating photon: loop limit reached."); @@ -57,20 +57,24 @@ pub fn raman( } // Interaction distances. - let voxel_dist = voxel - .dist(phot.ray()) - .expect("Could not determine voxel distance."); + let index = input.grid.gen_index_voxel(phot.ray().pos()); + let voxel_dist = match &index { + Some((_index, voxel)) => { + voxel.dist(phot.ray()).expect("Could not determine voxel distance.") + }, + None => f64::INFINITY, + }; let scat_dist = -(rng.gen::()).ln() / env.inter_coeff(); let surf_hit = input .tree .scan(phot.ray().clone(), bump_dist, voxel_dist.min(scat_dist)); - let boundary_hit = None; + let boundary_hit = input.bound.dist_boundary(phot.ray()).expect("Photon not contained in boundary. "); // Event handling. match Event::new(voxel_dist, scat_dist, surf_hit, boundary_hit, bump_dist) { - Event::Voxel(dist) => travel(&mut data, &mut phot, &env, index, dist + bump_dist), + Event::Voxel(dist) => travel(&mut data, &mut phot, &env, &index, dist + bump_dist), Event::Scattering(dist) => { - travel(&mut data, &mut phot, &env, index, dist); + travel(&mut data, &mut phot, &env, &index, dist); // // Capture. // if let Some(weight) = @@ -82,12 +86,12 @@ pub fn raman( shift_scatter(&mut rng, &mut phot, &env); } Event::Surface(hit) => { - travel(&mut data, &mut phot, &env, index, hit.dist()); + travel(&mut data, &mut phot, &env, &index, hit.dist()); surface(&mut rng, &hit, &mut phot, &mut env, &mut data); - travel(&mut data, &mut phot, &env, index, bump_dist); + travel(&mut data, &mut phot, &env, &index, bump_dist); }, Event::Boundary(boundary_hit) => { - travel(&mut data, &mut phot, &env, index, boundary_hit.dist()); + travel(&mut data, &mut phot, &env, &index, boundary_hit.dist()); input.bound.apply(rng, &boundary_hit, &mut phot); } } diff --git a/src/sim/engine/engines/standard.rs b/src/sim/engine/engines/standard.rs index bc93620..2aaee4b 100644 --- a/src/sim/engine/engines/standard.rs +++ b/src/sim/engine/engines/standard.rs @@ -1,7 +1,7 @@ //! Standard photon-lifetime engine function. use crate::{ - geom::{Boundary, Trace}, + geom::Trace, phys::Photon, sim::{scatter::scatter, surface::surface, travel::travel, Event, Input, Output}, }; @@ -31,7 +31,10 @@ pub fn standard(input: &Input, mut data: &mut Output, mut rng: &mut ThreadRng, m // Main event loop. let mut num_loops = 0; - while let Some((index, voxel)) = input.grid.gen_index_voxel(phot.ray().pos()) { + // It could be that this is preventing our photon packets from interacting with the boundary. + //while let Some((index, voxel)) = input.grid.gen_index_voxel(phot.ray().pos()) { + while input.bound.contains(phot.ray().pos()) { + // Loop limit check. if num_loops >= loop_limit { println!("[WARN] : Terminating photon: loop limit reached."); @@ -49,30 +52,35 @@ pub fn standard(input: &Input, mut data: &mut Output, mut rng: &mut ThreadRng, m } // Interaction distances. - let voxel_dist = voxel - .dist(phot.ray()) - .expect("Could not determine voxel distance."); + let index = input.grid.gen_index_voxel(phot.ray().pos()); + let voxel_dist = match &index { + Some((_index, voxel)) => { + voxel.dist(phot.ray()).expect("Could not determine voxel distance.") + }, + None => f64::INFINITY, + }; let scat_dist = -(rng.gen::()).ln() / env.inter_coeff(); let surf_hit = input .tree .scan(phot.ray().clone(), bump_dist, voxel_dist.min(scat_dist)); - let boundary_hit = None; + let boundary_hit = input.bound.dist_boundary(phot.ray()).expect("Photon not contained in boundary. "); // Event handling. match Event::new(voxel_dist, scat_dist, surf_hit, boundary_hit, bump_dist) { - Event::Voxel(dist) => travel(&mut data, &mut phot, &env, index, dist + bump_dist), + Event::Voxel(dist) => travel(&mut data, &mut phot, &env, &index, dist + bump_dist), Event::Scattering(dist) => { - travel(&mut data, &mut phot, &env, index, dist); + travel(&mut data, &mut phot, &env, &index, dist); scatter(&mut rng, &mut phot, &env); } Event::Surface(hit) => { - travel(&mut data, &mut phot, &env, index, hit.dist()); + travel(&mut data, &mut phot, &env, &index, hit.dist()); surface(&mut rng, &hit, &mut phot, &mut env, &mut data); - travel(&mut data, &mut phot, &env, index, bump_dist); + travel(&mut data, &mut phot, &env, &index, bump_dist); }, Event::Boundary(boundary_hit) => { - travel(&mut data, &mut phot, &env, index, boundary_hit.dist()); + travel(&mut data, &mut phot, &env, &index, boundary_hit.dist()); input.bound.apply(rng, &boundary_hit, &mut phot); + travel(&mut data, &mut phot, &env, &index, 100.0 * bump_dist); } } diff --git a/src/sim/event.rs b/src/sim/event.rs index 619d4d4..cb0813e 100644 --- a/src/sim/event.rs +++ b/src/sim/event.rs @@ -1,6 +1,6 @@ //! Event enumeration. -use crate::geom::{boundary, BoundaryHit, Hit}; +use crate::geom::{BoundaryHit, Hit}; /// Event determination enumeration. #[derive(PartialEq, Debug)] @@ -25,13 +25,15 @@ impl<'a, T> Event<'a, T> { voxel_dist: f64, scat_dist: f64, surf_hit: Option>, - boundary_hit: Option>, + boundary_hit: BoundaryHit<'a>, bump_dist: f64, ) -> Self { debug_assert!(voxel_dist > 0.0); debug_assert!(scat_dist > 0.0); debug_assert!(bump_dist > 0.0); + // Logically, if there is any geometry, it should be within the octree + // which is contained within the boundary. if let Some(hit) = surf_hit { if voxel_dist < (hit.dist() + bump_dist) { if scat_dist < (voxel_dist + bump_dist) { @@ -46,16 +48,10 @@ impl<'a, T> Event<'a, T> { return Self::Surface(hit); } - // We should be able to safely assume that if there were geometry in the - if let Some(bhit) = boundary_hit { - if bhit.dist() < scat_dist && bhit.dist() < (voxel_dist + bump_dist) { - return Self::Boundary(bhit) - } + if boundary_hit.dist() <= (voxel_dist + bump_dist) { + return Self::Boundary(boundary_hit); } - if scat_dist < (voxel_dist + bump_dist) { - return Self::Scattering(scat_dist); - } Self::Voxel(voxel_dist) } } @@ -64,7 +60,7 @@ impl<'a, T> Event<'a, T> { mod tests { use crate::{ sim::Attribute, - geom::{Side, BoundaryCondition}, + geom::{Side, BoundaryCondition, BoundaryDirection}, math::Dir3, }; use super::*; @@ -73,7 +69,8 @@ mod tests { #[test] fn test_new_surface_hit() { let surf_hit = Some(Hit::new(&Attribute::Mirror(0.5), 1.0, Side::Outside(Dir3::new(1.0, 0.0, 0.0)))); - let event = Event::new(2.0, 3.0, surf_hit, None, 0.5); + let boundary_hit = BoundaryHit::new(&BoundaryCondition::Kill, f64::INFINITY, BoundaryDirection::Bottom); + let event = Event::new(2.0, 3.0, surf_hit, boundary_hit, 0.5); // Check each of the components of the event. if let Event::Surface(hit) = event { @@ -87,22 +84,24 @@ mod tests { #[test] fn test_new_voxel_collision() { - let event: Event<'_, Attribute> = Event::new(2.0, 3.0, None, None, 0.5); + let boundary_hit = BoundaryHit::new(&BoundaryCondition::Kill, f64::INFINITY, BoundaryDirection::Bottom); + let event: Event<'_, Attribute> = Event::new(2.0, 3.0, None, boundary_hit, 0.5); assert_eq!(event, Event::Voxel(2.0)); } #[test] fn test_new_scattering_event() { let surf_hit = Some(Hit::new(&Attribute::Mirror(0.5), 2.0, Side::Outside(Dir3::new(1.0, 0.0, 0.0)))); - let event = Event::new(2.0, 1.0, surf_hit, None, 0.5); + let boundary_hit = BoundaryHit::new(&BoundaryCondition::Kill, f64::INFINITY, BoundaryDirection::Bottom); + let event = Event::new(2.0, 1.0, surf_hit, boundary_hit, 0.5); assert_eq!(event, Event::Scattering(1.0)); } #[test] fn test_new_boundary_event() { - let bhit = BoundaryHit::new(&BoundaryCondition::Periodic(0.0), 0.1, boundary::BoundaryDirection::North); - let event: Event<'_, Attribute<'_>> = Event::new(2.0, 1.0, None, Some(bhit.clone()), 0.5); + let bhit = BoundaryHit::new(&BoundaryCondition::Periodic(0.0), 0.1, BoundaryDirection::North); + let event: Event<'_, Attribute<'_>> = Event::new(2.0, 1.0, None, bhit.clone(), 0.5); assert_eq!(event, Event::Boundary(bhit)); } } diff --git a/src/sim/mod.rs b/src/sim/mod.rs index 94cc6e4..2b8b007 100644 --- a/src/sim/mod.rs +++ b/src/sim/mod.rs @@ -9,6 +9,7 @@ pub mod film_builder; pub mod frame; pub mod input; pub mod output; +pub mod output_volume; pub mod param; pub mod peel_off; pub mod photon_collector; @@ -19,6 +20,6 @@ pub mod surface; pub mod travel; pub use self::{ - attribute::*, engine::*, event::*, film_builder::*, frame::*, input::*, output::*, param::*, + attribute::*, engine::*, event::*, film_builder::*, frame::*, input::*, output::*, output_volume::*, param::*, peel_off::*, photon_collector::*, run::*, scatter::*, settings::*, surface::*, travel::*, }; diff --git a/src/sim/travel.rs b/src/sim/travel.rs index db6bbfa..ffc45f9 100644 --- a/src/sim/travel.rs +++ b/src/sim/travel.rs @@ -3,18 +3,24 @@ use crate::{ phys::{Local, Photon}, sim::Output, + geom::Cube, }; use physical_constants::SPEED_OF_LIGHT_IN_VACUUM; /// Move the photon forward and record the flight. #[inline] -pub fn travel(data: &mut Output, phot: &mut Photon, env: &Local, index: [usize; 3], dist: f64) { +pub fn travel(data: &mut Output, phot: &mut Photon, env: &Local, index: &Option<([usize; 3], Cube)>, dist: f64) { debug_assert!(dist > 0.0); - let weight_power_dist = phot.weight() * phot.power() * dist; - data.energy[index] += weight_power_dist * env.ref_index() / SPEED_OF_LIGHT_IN_VACUUM; - data.absorptions[index] += weight_power_dist * env.abs_coeff(); - data.shifts[index] += weight_power_dist * env.shift_coeff(); + match *index { + Some((idx, _)) => { + let weight_power_dist = phot.weight() * phot.power() * dist; + data.energy[idx] += weight_power_dist * env.ref_index() / SPEED_OF_LIGHT_IN_VACUUM; + data.absorptions[idx] += weight_power_dist * env.abs_coeff(); + data.shifts[idx] += weight_power_dist * env.shift_coeff(); + }, + None => {}, + } phot.ray_mut().travel(dist); } From fa9c6cd600bc21d328e86f4f1e9247c8e3efe60a Mon Sep 17 00:00:00 2001 From: Sam Morrell Date: Thu, 18 Apr 2024 09:37:10 +0100 Subject: [PATCH 03/16] Removed import of not-yet-existing module. --- src/sim/mod.rs | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/src/sim/mod.rs b/src/sim/mod.rs index 2b8b007..94cc6e4 100644 --- a/src/sim/mod.rs +++ b/src/sim/mod.rs @@ -9,7 +9,6 @@ pub mod film_builder; pub mod frame; pub mod input; pub mod output; -pub mod output_volume; pub mod param; pub mod peel_off; pub mod photon_collector; @@ -20,6 +19,6 @@ pub mod surface; pub mod travel; pub use self::{ - attribute::*, engine::*, event::*, film_builder::*, frame::*, input::*, output::*, output_volume::*, param::*, + attribute::*, engine::*, event::*, film_builder::*, frame::*, input::*, output::*, param::*, peel_off::*, photon_collector::*, run::*, scatter::*, settings::*, surface::*, travel::*, }; From 59c4f43ab36666c64d2eda2ad01ed217b2f07ea4 Mon Sep 17 00:00:00 2001 From: Sam Morrell Date: Thu, 18 Apr 2024 19:22:56 +0100 Subject: [PATCH 04/16] Fixed 'coneditions' typo in comment. --- src/bin/mcrt.rs | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/bin/mcrt.rs b/src/bin/mcrt.rs index 72d92bb..f66c098 100644 --- a/src/bin/mcrt.rs +++ b/src/bin/mcrt.rs @@ -85,7 +85,7 @@ fn main() { report!(surfs, "surfaces"); /* - * Create a boundary for the simulation with boundary coneditions. + * Create a boundary for the simulation with boundary conditions. * For now we hard-code this to kill, but we can link this to configuration soon. * TODO: We probably want to implement the MPI adjacent rank transfer here too. */ From 5edde8ef896711f4f528066d178014ddf942ae40 Mon Sep 17 00:00:00 2001 From: Sam Morrell Date: Thu, 1 Aug 2024 16:05:10 +0100 Subject: [PATCH 05/16] Added new output types - Started refactor of output with new `io` module. - Added volume output + tests. - Added plane output + tests. --- Cargo.toml | 4 + src/io/mod.rs | 5 + src/io/output/mod.rs | 5 + src/io/output/output_plane.rs | 186 +++++++++++++++++++++++++++++++++ src/io/output/output_volume.rs | 79 ++++++++++++++ src/lib.rs | 1 + src/sim/output.rs | 2 +- 7 files changed, 281 insertions(+), 1 deletion(-) create mode 100644 src/io/mod.rs create mode 100644 src/io/output/mod.rs create mode 100644 src/io/output/output_plane.rs create mode 100644 src/io/output/output_volume.rs diff --git a/Cargo.toml b/Cargo.toml index 1e394ee..61cc9d9 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -15,6 +15,7 @@ exclude = ["docs/", "input/", "output/"] [features] default = [] mpi = ["dep:mpi", "dep:memoffset"] +pyo3 = ['dep:pyo3'] [dependencies] # ARCTK Dependencies @@ -56,6 +57,9 @@ netcdf = "^0.8.1" mpi = { version = "0.7.0", optional = true } memoffset = { version = "0.9.0", optional = true } +# Python modules +pyo3 = { version = "0.21.1", features = ["extension-module"], optional = true } + [dev-dependencies] tempfile = "3.2.*" statrs = "0.15.*" diff --git a/src/io/mod.rs b/src/io/mod.rs new file mode 100644 index 0000000..5f52c48 --- /dev/null +++ b/src/io/mod.rs @@ -0,0 +1,5 @@ +/** + * The I/O module. + * This contains for code the input and outputs of an MCRT simulation. + */ +pub mod output; \ No newline at end of file diff --git a/src/io/output/mod.rs b/src/io/output/mod.rs new file mode 100644 index 0000000..a84e455 --- /dev/null +++ b/src/io/output/mod.rs @@ -0,0 +1,5 @@ + +pub mod output_plane; +pub mod output_volume; + +pub use self::{output_plane::*, output_volume::*}; \ No newline at end of file diff --git a/src/io/output/output_plane.rs b/src/io/output/output_plane.rs new file mode 100644 index 0000000..e131796 --- /dev/null +++ b/src/io/output/output_plane.rs @@ -0,0 +1,186 @@ +use ndarray::Array2; +use crate::{ + math::Point2, + ord::cartesian::{X, Y} +}; + +pub enum AxisAlignedPlane{ + XY, + XZ, + YZ, +} + +#[derive(Debug)] +pub struct OutputPlane { + mins: Point2, + maxs: Point2, + res: [usize; 2], + data: Array2, +} + +impl OutputPlane { + pub fn new(mins: Point2, maxs: Point2, res: [usize; 2]) -> Self { + debug_assert!(res[X] > 0); + debug_assert!(res[Y] > 0); + + Self { + mins, + maxs, + res, + data: Array2::zeros(res) + } + } + + pub fn xlen(&self) -> f64 { + self.maxs.x() - self.mins.x() + } + + pub fn ylen(&self) -> f64 { + self.maxs.y() - self.mins.y() + } + + pub fn dx(&self) -> f64 { + self.xlen() / self.res[X] as f64 + } + + pub fn dy(&self) -> f64 { + self.ylen() / self.res[Y] as f64 + } + + pub fn area(&self) -> f64 { + (self.maxs.x() - self.mins.x()) * (self.maxs.y() - self.mins.y()) + } + + pub fn pix_area(&self) -> f64 { + self.area() / (self.res[X] * self.res[Y]) as f64 + } + + pub fn index_for_coord(&self, x: f64, y:f64) -> Option<(usize, usize)> { + if x < self.mins[X] { return None; } + let xoff = x - self.mins[X] as f64; + if y < self.mins[Y] { return None; } + let yoff = y - self.mins[Y] as f64; + + if xoff < self.xlen() && yoff < self.ylen() { + let i = (xoff / self.dx()).floor() as usize; + let j = (yoff / self.dy()).floor() as usize; + Some((i, j)) + } else { + None + } + } + + pub fn at(&self, x: f64, y: f64) -> Option<&f64> { + let (i, j) = self.index_for_coord(x, y)?; + Some(&self.data[[i, j]]) + } + + pub fn at_mut(&mut self, x: f64, y: f64) -> Option<&mut f64> { + let (i, j) = self.index_for_coord(x, y)?; + Some(&mut self.data[[i, j]]) + } + +} + +#[cfg(test)] +mod tests { + use super::*; + + #[test] + fn test_data_allocation() { + let mins = Point2::new(0.0, 0.0); + let maxs = Point2::new(10.0, 10.0); + let res = [100, 100]; + let output_plane = OutputPlane::new(mins, maxs, res); + + // Check that the data array is correctly allocated with the specified resolution + assert_eq!(output_plane.data.shape(), &[100, 100]); + assert_eq!(output_plane.data.len(), 10000); + } + + // TODO: Check this test. + #[test] + fn test_area() { + let mins = Point2::new(0.0, 0.0); + let maxs = Point2::new(10.0, 10.0); + let res = [100, 100]; + let output_plane = OutputPlane::new(mins, maxs, res); + + let expected_area = 100.0; + let actual_area = output_plane.area(); + + assert_eq!(actual_area, expected_area); + } + + #[test] + fn test_dx_dy() { + let mins = Point2::new(0.0, 0.0); + let maxs = Point2::new(10.0, 5.0); + let res = [100, 100]; + let output_plane = OutputPlane::new(mins, maxs, res); + + // Check the result based on the dimensions. + assert_eq!(output_plane.dx(), 0.1); + assert_eq!(output_plane.dy(), 0.05); + } + + // TODO: Check this test. + #[test] + fn test_pix_area() { + let mins = Point2::new(0.0, 0.0); + let maxs = Point2::new(10.0, 10.0); + let res = [100, 100]; + let output_plane = OutputPlane::new( mins, maxs, res); + let expected_pix_area = 0.01; + let actual_pix_area = output_plane.pix_area(); + assert_eq!(actual_pix_area, expected_pix_area); + } + + #[test] + fn test_index_for_coord() { + let mins = Point2::new(0.0, 0.0); + let maxs = Point2::new(10.0, 10.0); + let res = [100, 100]; + let output_plane = OutputPlane::new(mins, maxs, res); + + // Check the pixel indices + assert_eq!(output_plane.index_for_coord(5.0, 5.0), Some((50, 50))); + assert_eq!(output_plane.index_for_coord(7.0, 7.0), Some((70, 70))); + assert_eq!(output_plane.index_for_coord(9.0, 9.0), Some((90, 90))); + assert_eq!(output_plane.index_for_coord(0.0, 0.0), Some((0, 0))); + assert_eq!(output_plane.index_for_coord(10.0, 10.0), None); // Outside the grid, outer extremity. + } + + #[test] + fn test_index_for_coord_outside_grid() { + let mins = Point2::new(0.0, 0.0); + let maxs = Point2::new(10.0, 10.0); + let res = [100, 100]; + let output_plane = OutputPlane::new(mins, maxs, res); + + // Test coordinates outside of the grid + assert_eq!(output_plane.index_for_coord(-1.0, 5.0), None); + assert_eq!(output_plane.index_for_coord(5.0, -1.0), None); + assert_eq!(output_plane.index_for_coord(11.0, 5.0), None); + assert_eq!(output_plane.index_for_coord(5.0, 11.0), None); + } + + #[test] + #[ignore] + fn test_edit_value() { + let mins = Point2::new(0.0, 0.0); + let maxs = Point2::new(10.0, 10.0); + let res = [100, 100]; + let mut output_plane = OutputPlane::new(mins, maxs, res); + + // Set initial values + *output_plane.at_mut(5.0, 5.0).unwrap() = 1.0; + *output_plane.at_mut(7.0, 7.0).unwrap() = 2.0; + *output_plane.at_mut(9.0, 9.0).unwrap() = 3.0; + + // Check that the values have been updated + assert_eq!(output_plane.at(5.0, 5.0), Some(&1.0)); + assert_eq!(output_plane.at(7.0, 7.0), Some(&2.0)); + assert_eq!(output_plane.at(9.0, 9.0), Some(&3.0)); + } +} \ No newline at end of file diff --git a/src/io/output/output_volume.rs b/src/io/output/output_volume.rs new file mode 100644 index 0000000..aae8f4e --- /dev/null +++ b/src/io/output/output_volume.rs @@ -0,0 +1,79 @@ +use std::ops::AddAssign; +#[cfg(feature = "pyo3")] +use pyo3::prelude::*; + +use ndarray::Array3; +use crate::{ + access, fs::Save, geom::Cube, ord::cartesian::{X, Y, Z} +}; + +#[derive(PartialEq, Debug, Clone)] +#[cfg_attr(feature = "pyo3", pyclass)] +pub enum OutputParameter { + Energy, + Absorption, + Shift, +} + +#[derive(Debug)] +pub struct OutputVolume { + boundary: Cube, + res: [usize; 3], + param: OutputParameter, + data: Array3, +} + +impl OutputVolume { + access!(boundary: Cube); + access!(res: [usize; 3]); + access!(data, data_mut: Array3); + + pub fn new(boundary: Cube, res: [usize; 3], param: OutputParameter) -> Self { + // Check that we don't have non-zero cell size. + debug_assert!(res[X] > 0); + debug_assert!(res[Y] > 0); + debug_assert!(res[Z] > 0); + + Self { + boundary, + res, + param, + data: Array3::zeros(res) + } + } + + pub fn cell_volume(&self) -> f64 { + self.boundary.vol() / (self.res[X] * self.res[Y] * self.res[Z]) as f64 + } +} + +impl AddAssign<&Self> for OutputVolume { + fn add_assign(&mut self, rhs: &Self) { + debug_assert_eq!(self.boundary(), rhs.boundary()); + debug_assert_eq!(self.cell_volume(), rhs.cell_volume()); + debug_assert_eq!(self.param, rhs.param); + + self.data += &rhs.data + } +} + +impl Save for OutputVolume { + #[inline] + fn save_data(&self, path: &std::path::Path) -> Result<(), crate::err::Error> { + (&self.data / self.cell_volume()).save(&path)?; + Ok(()) + } +} + +#[cfg(test)] +mod tests { + use crate::math::Point3; + use super::*; + + #[test] + fn new_test() { + let boundary = Cube::new(Point3::new(0.0,0.0,0.0), Point3::new(10.0,10.0,10.0)); + let outvol = OutputVolume::new(boundary, [10, 10, 10], OutputParameter::Energy); + assert_eq!(outvol.cell_volume(), 1.0); + } +} \ No newline at end of file diff --git a/src/lib.rs b/src/lib.rs index f050d26..8e27495 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -41,6 +41,7 @@ pub mod err; pub mod fs; pub mod geom; pub mod img; +pub mod io; pub mod math; pub mod meta; pub mod ord; diff --git a/src/sim/output.rs b/src/sim/output.rs index 6a01360..e144915 100644 --- a/src/sim/output.rs +++ b/src/sim/output.rs @@ -33,7 +33,7 @@ pub struct Output<'a> { /// Emission power. pub emission: Array3, - /// Photo-energy. + /// Photon-energy. pub energy: Array3, /// Absorptions. pub absorptions: Array3, From a1f055896208948069edf7a8424c35583875a55b Mon Sep 17 00:00:00 2001 From: Sam Morrell Date: Tue, 13 Aug 2024 09:29:23 +0100 Subject: [PATCH 06/16] Moved Photon collector to `io::output`. It makes more sense for the `PhotonCollector` struct to be contained in the output module, so I have moved it there. I have changed references in the code to point to this new location. --- src/bin/mcrt.rs | 3 +- src/io/output/mod.rs | 7 +++- src/{sim => io/output}/photon_collector.rs | 38 ----------------- src/io/output/photon_collector_builder.rs | 49 ++++++++++++++++++++++ src/sim/mod.rs | 3 +- src/sim/output.rs | 2 +- 6 files changed, 59 insertions(+), 43 deletions(-) rename src/{sim => io/output}/photon_collector.rs (78%) create mode 100644 src/io/output/photon_collector_builder.rs diff --git a/src/bin/mcrt.rs b/src/bin/mcrt.rs index f66c098..39feae2 100644 --- a/src/bin/mcrt.rs +++ b/src/bin/mcrt.rs @@ -14,9 +14,10 @@ use aetherus::{ img::{Colour, Image}, ord::{Build, Link, Register, Set, cartesian::{X, Y}}, report, + io::output::PhotonCollector, sim::{ run, AttributeLinkerLinkerLinkerLinkerLinker as Attr, Engine, Input, Output, Parameters, - ParametersBuilderLoader, PhotonCollector, + ParametersBuilderLoader, }, util::{ banner::{section, sub_section, title}, diff --git a/src/io/output/mod.rs b/src/io/output/mod.rs index a84e455..22cd20d 100644 --- a/src/io/output/mod.rs +++ b/src/io/output/mod.rs @@ -1,5 +1,10 @@ pub mod output_plane; pub mod output_volume; +pub mod photon_collector; +pub mod photon_collector_builder; -pub use self::{output_plane::*, output_volume::*}; \ No newline at end of file +pub use self::{ + photon_collector::*, + photon_collector_builder::*, +}; \ No newline at end of file diff --git a/src/sim/photon_collector.rs b/src/io/output/photon_collector.rs similarity index 78% rename from src/sim/photon_collector.rs rename to src/io/output/photon_collector.rs index 799904c..fb94962 100644 --- a/src/sim/photon_collector.rs +++ b/src/io/output/photon_collector.rs @@ -1,44 +1,6 @@ use crate::{err::Error, fmt_report, fs::Save, phys::Photon, tools::ProgressBar}; use std::{fmt::Display, fs::File, io::Write, ops::AddAssign, path::Path}; -/* -pub fn photon_collector_write_thread(path_str: String, rx: Receiver>) { - let file_path = Path::new(&path_str); - let mut file = File::create(file_path).expect("Unable to open file for writing. "); - - // To reduce the time to run, I am manually do my won CSV write, directly from this vec. - let headings = vec!["pos_x", "pos_y", "pos_z", "dir_x", "dir_y", "dir_z", "wavelength", "weight"]; - write!(file, "{}", headings[0]).expect("Unable to write header to file. ");; - for heading in headings.iter().skip(1) { - write!(file, ",{}", heading).expect("Unable to write header to file. "); - } - writeln!(file).expect("Unable to write header to file. "); - - loop { - while let optphot = rx.recv().unwrap() { - match optphot { - Some(phot) => { - writeln!(file, "{},{},{},{},{},{},{},{}", - phot.ray().pos().x(), - phot.ray().pos().y(), - phot.ray().pos().z(), - phot.ray().dir().x(), - phot.ray().dir().y(), - phot.ray().dir().z(), - phot.wavelength(), - phot.weight(), - ).expect("Unable to write photon to file. "); - } - None => break, - } - } - - // Wait for new photons. - thread::sleep(Duration::from_millis(1000)); - } -} -*/ - #[derive(Default, Clone)] pub struct PhotonCollector { /// The vector of collected photons. diff --git a/src/io/output/photon_collector_builder.rs b/src/io/output/photon_collector_builder.rs new file mode 100644 index 0000000..b5dc20a --- /dev/null +++ b/src/io/output/photon_collector_builder.rs @@ -0,0 +1,49 @@ +use serde::{Serialize, Deserialize}; +use crate::io::output::PhotonCollector; + +#[derive(Default, Debug, Serialize, Deserialize)] +pub struct PhotonCollectorBuilder { + kill_photons: Option, +} + +impl PhotonCollectorBuilder { + pub fn build(&self) -> PhotonCollector { + let mut photcol = PhotonCollector::new(); + photcol.kill_photon = match self.kill_photons { + Some(kp) => kp, + None => false + }; + photcol + } +} + +#[cfg(test)] +mod tests { + use super::*; + use json5; + + #[test] + fn test_new() { + let photcol = PhotonCollectorBuilder::default().build(); + assert_eq!(photcol.kill_photon, false); + assert_eq!(photcol.nphoton(), 0); + } + + #[test] + fn test_deserialise_default() { + let input = "{}"; + let photcolbuild: PhotonCollectorBuilder = json5::from_str(&input).unwrap(); + let photcol = photcolbuild.build(); + assert_eq!(photcol.kill_photon, false); + assert_eq!(photcol.nphoton(), 0); + } + + #[test] + fn test_deserialise_kill_photons() { + let input = "{ kill_photons: true }"; + let photcolbuild: PhotonCollectorBuilder = json5::from_str(&input).unwrap(); + let photcol = photcolbuild.build(); + assert_eq!(photcol.kill_photon, true); + assert_eq!(photcol.nphoton(), 0); + } +} \ No newline at end of file diff --git a/src/sim/mod.rs b/src/sim/mod.rs index 94cc6e4..083b4b3 100644 --- a/src/sim/mod.rs +++ b/src/sim/mod.rs @@ -11,7 +11,6 @@ pub mod input; pub mod output; pub mod param; pub mod peel_off; -pub mod photon_collector; pub mod run; pub mod scatter; pub mod settings; @@ -20,5 +19,5 @@ pub mod travel; pub use self::{ attribute::*, engine::*, event::*, film_builder::*, frame::*, input::*, output::*, param::*, - peel_off::*, photon_collector::*, run::*, scatter::*, settings::*, surface::*, travel::*, + peel_off::*, run::*, scatter::*, settings::*, surface::*, travel::*, }; diff --git a/src/sim/output.rs b/src/sim/output.rs index e144915..e6a8306 100644 --- a/src/sim/output.rs +++ b/src/sim/output.rs @@ -13,6 +13,7 @@ use crate::{ cartesian::{X, Y, Z} }, util::fmt::DataCube, + io::output::PhotonCollector, }; use ndarray::Array3; use std::{ @@ -21,7 +22,6 @@ use std::{ path::Path, }; -use super::PhotonCollector; /// MCRT output data. #[derive(Clone)] From dbfc01a343392715ed8bd38f8126f88b4ca16e40 Mon Sep 17 00:00:00 2001 From: Sam Morrell Date: Tue, 13 Aug 2024 09:32:47 +0100 Subject: [PATCH 07/16] Added `Output` struct in `io::output`. This struct is intended as the new central object containing all of the output code for `mcrt`. It contains all of the other outputs. I am in the process of refactoring all of the output structs to work with this. - Added `Output` struct. - Added builder structs which support deserialisation from JSON to allow init from config. - Added supporting structs around output types to deserialise / integrate. --- src/io/output/mod.rs | 15 +++- src/io/output/output.rs | 27 +++++++ src/io/output/output_builder.rs | 5 ++ src/io/output/output_config.rs | 107 +++++++++++++++++++++++++ src/io/output/output_plane.rs | 18 ++++- src/io/output/output_plane_builder.rs | 38 +++++++++ src/io/output/output_registry.rs | 17 ++++ src/io/output/output_type.rs | 9 +++ src/io/output/output_volume.rs | 9 ++- src/io/output/output_volume_builder.rs | 43 ++++++++++ 10 files changed, 282 insertions(+), 6 deletions(-) create mode 100644 src/io/output/output.rs create mode 100644 src/io/output/output_builder.rs create mode 100644 src/io/output/output_config.rs create mode 100644 src/io/output/output_plane_builder.rs create mode 100644 src/io/output/output_registry.rs create mode 100644 src/io/output/output_type.rs create mode 100644 src/io/output/output_volume_builder.rs diff --git a/src/io/output/mod.rs b/src/io/output/mod.rs index 22cd20d..63ef456 100644 --- a/src/io/output/mod.rs +++ b/src/io/output/mod.rs @@ -1,10 +1,23 @@ - +pub mod output; +pub mod output_config; pub mod output_plane; +pub mod output_plane_builder; +pub mod output_registry; +pub mod output_type; pub mod output_volume; +pub mod output_volume_builder; pub mod photon_collector; pub mod photon_collector_builder; pub use self::{ + output::*, + output_config::*, + output_plane::*, + output_plane_builder::*, + output_registry::*, + output_type::*, + output_volume::*, + output_volume_builder::*, photon_collector::*, photon_collector_builder::*, }; \ No newline at end of file diff --git a/src/io/output/output.rs b/src/io/output/output.rs new file mode 100644 index 0000000..f87fff1 --- /dev/null +++ b/src/io/output/output.rs @@ -0,0 +1,27 @@ +use crate::{ + data::Histogram, + img::Image, + io::output::{OutputPlane, OutputVolume, OutputRegistry, PhotonCollector}, +}; +use ndarray::Array3; + +pub struct Output { + /// Output volumes. + pub vol: Vec, + /// Output planes. + pub plane: Vec, + /// Photon Collectors. + pub phot_cols: Vec, + /// Spectra. + pub specs: Vec, + /// Image data. + pub imgs: Vec, + /// CCD Data. + pub ccds: Vec>, + /// Photo data. + pub photos: Vec, + + /// Contains the mapping between index and name for + /// each of the output types. + pub reg: OutputRegistry, +} \ No newline at end of file diff --git a/src/io/output/output_builder.rs b/src/io/output/output_builder.rs new file mode 100644 index 0000000..2a724c9 --- /dev/null +++ b/src/io/output/output_builder.rs @@ -0,0 +1,5 @@ +use crate::io::Output; + +pub struct OutputBuilder { + +} \ No newline at end of file diff --git a/src/io/output/output_config.rs b/src/io/output/output_config.rs new file mode 100644 index 0000000..404f87c --- /dev/null +++ b/src/io/output/output_config.rs @@ -0,0 +1,107 @@ +use std::{ + collections::HashMap, + fmt +}; + +use serde::{Serialize}; +use arctk_attr::file; +use crate::{io::output::{OutputPlaneBuilder, OutputVolumeBuilder, PhotonCollectorBuilder}}; + + +#[file] +#[derive(Serialize)] +pub struct OutputConfig { + pub volumes: Option>, + pub planes: Option>, + pub photon_collectors: Option>, +} + +impl OutputConfig { + pub fn n_volumes(&self) -> usize { + match &self.volumes { + Some(vol) => vol.iter().count(), + None => 0, + } + } + + pub fn n_planes(&self) -> usize { + match &self.planes { + Some(plane) => plane.iter().count(), + None => 0, + } + } + + pub fn n_photon_collectors(&self) -> usize { + match &self.photon_collectors { + Some(pc) => pc.iter().count(), + None => 0, + } + } +} + +impl fmt::Display for OutputConfig { + fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result { + write!(f, "OutputConfig {{\n")?; + + if let Some(volumes) = &self.volumes { + write!(f, " volumes: {{\n")?; + for (name, volume) in volumes { + write!(f, " {}: {:?}\n", name, volume)?; + } + write!(f, " }}\n")?; + } + + if let Some(planes) = &self.planes { + write!(f, " planes: {{\n")?; + for (name, plane) in planes { + write!(f, " {}: {:?}\n", name, plane)?; + } + write!(f, " }}\n")?; + } + + if let Some(photon_collectors) = &self.photon_collectors { + write!(f, " photon_collectors: {{\n")?; + for (name, collector) in photon_collectors { + write!(f, " {}: {:?}\n", name, collector)?; + } + write!(f, " }}\n")?; + } + + write!(f, "}}") + } +} + +#[cfg(test)] +mod tests { + + use super::*; + use json5; + + #[test] + fn output_config_deserialise_test() { + let conf_str = r#" + { + volumes: { + full_vol: { boundary: [[0, 0, 0], [10, 10, 10]], res: [10, 10, 10], param: "energy" }, + partial_vol: { boundary: [[2.5, 2.5, 0], [2.5, 2.5, 10]], res: [100, 100, 10], param: "energy" }, + }, + planes: { + bottom: { boundary: [[0, 0], [10, 10]], res: [10, 10], plane: "xy" }, + }, + photon_collectors: { + terrain_collector: { kill_photons: false }, + sky_collector: { kill_photons: true }, + } + } + "#; + + // Deserialise from the provided string above. + let conf: OutputConfig = json5::from_str(conf_str).unwrap(); + + // Check that all outputs make it through. + assert_eq!(conf.n_volumes(), 2); + assert_eq!(conf.n_planes(), 1); + assert_eq!(conf.n_photon_collectors(), 2); + } + +} \ No newline at end of file diff --git a/src/io/output/output_plane.rs b/src/io/output/output_plane.rs index e131796..06b59da 100644 --- a/src/io/output/output_plane.rs +++ b/src/io/output/output_plane.rs @@ -1,15 +1,29 @@ use ndarray::Array2; +use serde::{Deserialize, Serialize}; use crate::{ - math::Point2, + math::{Point2, Point3}, ord::cartesian::{X, Y} }; -pub enum AxisAlignedPlane{ +#[derive(PartialEq, Debug, Clone, Serialize, Deserialize)] +#[cfg_attr(feature = "pyo3", pyclass)] +#[serde(rename_all = "lowercase")] +pub enum AxisAlignedPlane { XY, XZ, YZ, } +impl AxisAlignedPlane { + pub fn project_onto_plane(&self, p: &Point3) -> (f64, f64) { + match *self { + Self::XY => (p.x(), p.y()), + Self::XZ => (p.x(), p.z()), + Self::YZ => (p.y(), p.z()), + } + } +} + #[derive(Debug)] pub struct OutputPlane { mins: Point2, diff --git a/src/io/output/output_plane_builder.rs b/src/io/output/output_plane_builder.rs new file mode 100644 index 0000000..d2ef79a --- /dev/null +++ b/src/io/output/output_plane_builder.rs @@ -0,0 +1,38 @@ +use serde::{Deserialize, Serialize}; +use crate::{ + math::Point2, + io::output::{OutputPlane, AxisAlignedPlane}, +}; + +#[derive(Debug, Deserialize, Serialize)] +pub struct OutputPlaneBuilder { + boundary: (Point2, Point2), + res: [usize; 2], + plane: AxisAlignedPlane, +} + +impl OutputPlaneBuilder { + pub fn build(&self) -> OutputPlane { + OutputPlane::new(self.boundary.0, self.boundary.1, self.res) + } +} + +#[cfg(test)] +mod tests { + use super::*; + + #[test] + fn new_deserialise_build() { + let input = r#" + { + plane: "xy", + boundary: [[0, 0], [10, 10]], + res: [10, 10], + } + "#; + + let builder: OutputPlaneBuilder = json5::from_str(input).unwrap(); + let outvol = builder.build(); + assert_eq!(outvol.pix_area(), 1.0); + } +} diff --git a/src/io/output/output_registry.rs b/src/io/output/output_registry.rs new file mode 100644 index 0000000..c16efd9 --- /dev/null +++ b/src/io/output/output_registry.rs @@ -0,0 +1,17 @@ +use crate::{ + ord::Register, + io::output::Output, +}; + +pub struct OutputRegistry { + pub vol_reg: Register, + pub plane_reg: Register, + pub phot_cols_reg: Register, + pub spec_reg: Register, + pub img_reg: Register, + pub ccd_reg: Register, +} + +impl OutputRegistry { + +} \ No newline at end of file diff --git a/src/io/output/output_type.rs b/src/io/output/output_type.rs new file mode 100644 index 0000000..908aaed --- /dev/null +++ b/src/io/output/output_type.rs @@ -0,0 +1,9 @@ +pub enum OutputType { + Volume, + Plane, + PhotonCollector, + Spectrum, + Image, + Ccd, + Photo, +} \ No newline at end of file diff --git a/src/io/output/output_volume.rs b/src/io/output/output_volume.rs index aae8f4e..c81ed72 100644 --- a/src/io/output/output_volume.rs +++ b/src/io/output/output_volume.rs @@ -3,12 +3,14 @@ use std::ops::AddAssign; use pyo3::prelude::*; use ndarray::Array3; +use serde::{Deserialize, Serialize}; use crate::{ - access, fs::Save, geom::Cube, ord::cartesian::{X, Y, Z} + access, core::Real, fs::Save, geom::Cube, ord::cartesian::{X, Y, Z} }; -#[derive(PartialEq, Debug, Clone)] +#[derive(PartialEq, Debug, Clone, Serialize, Deserialize)] #[cfg_attr(feature = "pyo3", pyclass)] +#[serde(rename_all = "lowercase")] pub enum OutputParameter { Energy, Absorption, @@ -68,6 +70,7 @@ impl Save for OutputVolume { #[cfg(test)] mod tests { use crate::math::Point3; + use json5; use super::*; #[test] @@ -76,4 +79,4 @@ mod tests { let outvol = OutputVolume::new(boundary, [10, 10, 10], OutputParameter::Energy); assert_eq!(outvol.cell_volume(), 1.0); } -} \ No newline at end of file +} diff --git a/src/io/output/output_volume_builder.rs b/src/io/output/output_volume_builder.rs new file mode 100644 index 0000000..184b3fc --- /dev/null +++ b/src/io/output/output_volume_builder.rs @@ -0,0 +1,43 @@ +use serde::{Deserialize, Serialize}; +use crate::{ + geom::Cube, + math::Vec3, + io::output::{OutputVolume, OutputParameter}, +}; + +/// Configuration for the OutputVolume. +/// Importantly this can be serialised / deserialised using serde, so that this +/// can be built from a JSON configuration file. +#[derive(Serialize, Deserialize, Debug)] +pub struct OutputVolumeBuilder { + boundary: (Vec3, Vec3), + res: [usize; 3], + param: OutputParameter, +} + +impl OutputVolumeBuilder { + pub fn build(&self) -> OutputVolume { + let bound = Cube::new(self.boundary.0.data().into(), self.boundary.1.data().into()); + OutputVolume::new(bound, self.res, self.param.clone()) + } +} + +#[cfg(test)] +mod tests { + use super::*; + + #[test] + fn new_deserialise_build() { + let input = r#" + { + boundary: [[0, 0, 0], [10, 10, 10]], + res: [10, 10, 10], + param: "energy", + } + "#; + + let builder: OutputVolumeBuilder = json5::from_str(input).unwrap(); + let outvol = builder.build(); + assert_eq!(outvol.cell_volume(), 1.0); + } +} From 660c22a6b1df06b47dcf68c5be927cfd2d665e23 Mon Sep 17 00:00:00 2001 From: Sam Morrell Date: Wed, 14 Aug 2024 01:19:22 +0100 Subject: [PATCH 08/16] Added first implementation of new output and builder. - Output struct in `io` module now implemented. - Added Builder structure to deserialise and build Output. - Added supporting deserialisation, config and builder structs for output types. - Added voxel distance calculation to output volume to support next steps. - Added `OutputRegistry` to support output linking. --- src/data/histogram_builder.rs | 17 +++ src/data/mod.rs | 3 +- src/img/image_builder.rs | 14 +++ src/img/mod.rs | 3 +- src/io/output/ccd_builder.rs | 15 +++ src/io/output/mod.rs | 2 + src/io/output/output.rs | 19 ++- src/io/output/output_config.rs | 200 ++++++++++++++++++++++++++++++- src/io/output/output_registry.rs | 39 +++++- src/io/output/output_volume.rs | 66 +++++++++- 10 files changed, 364 insertions(+), 14 deletions(-) create mode 100644 src/data/histogram_builder.rs create mode 100644 src/img/image_builder.rs create mode 100644 src/io/output/ccd_builder.rs diff --git a/src/data/histogram_builder.rs b/src/data/histogram_builder.rs new file mode 100644 index 0000000..81abbb5 --- /dev/null +++ b/src/data/histogram_builder.rs @@ -0,0 +1,17 @@ +use serde::{Deserialize, Serialize}; + +use super::Histogram; + + +#[derive(Debug, Serialize, Deserialize)] +pub struct HistogramBuilder { + min: f64, + max: f64, + bins: usize, +} + +impl HistogramBuilder { + pub fn build(&self) -> Histogram { + Histogram::new(self.min, self.max, self.bins) + } +} \ No newline at end of file diff --git a/src/data/mod.rs b/src/data/mod.rs index 3bace4c..0fb70b5 100644 --- a/src/data/mod.rs +++ b/src/data/mod.rs @@ -2,7 +2,8 @@ pub mod average; pub mod histogram; +pub mod histogram_builder; pub mod histogram_iter; pub mod table; -pub use self::{average::*, histogram::*, histogram_iter::*, table::*}; +pub use self::{average::*, histogram::*, histogram_builder::*, histogram_iter::*, table::*}; diff --git a/src/img/image_builder.rs b/src/img/image_builder.rs new file mode 100644 index 0000000..df33038 --- /dev/null +++ b/src/img/image_builder.rs @@ -0,0 +1,14 @@ +use serde::{Serialize, Deserialize}; +use crate::img::{Image, Colour}; + +#[derive(Debug, Serialize, Deserialize)] +pub struct ImageBuilder { + res: [usize; 2], +} + +impl ImageBuilder { + pub fn build(&self) -> Image { + let background = Colour::new(0.0, 0.0, 0.0, 1.0); + Image::new_blank(self.res, background) + } +} \ No newline at end of file diff --git a/src/img/mod.rs b/src/img/mod.rs index c7e0446..1f1c9c0 100644 --- a/src/img/mod.rs +++ b/src/img/mod.rs @@ -5,5 +5,6 @@ pub mod colour; pub mod gradient; pub mod gradient_builder; pub mod image; +pub mod image_builder; -pub use self::{aspect_ratio::*, colour::*, gradient::*, gradient_builder::*, image::*}; +pub use self::{aspect_ratio::*, colour::*, gradient::*, gradient_builder::*, image::*, image_builder::*}; diff --git a/src/io/output/ccd_builder.rs b/src/io/output/ccd_builder.rs new file mode 100644 index 0000000..a378ba7 --- /dev/null +++ b/src/io/output/ccd_builder.rs @@ -0,0 +1,15 @@ +use serde::{Serialize, Deserialize}; +use ndarray::Array3; +use crate::ord::cartesian::{X, Y}; + +#[derive(Debug, Serialize, Deserialize)] +pub struct CcdBuilder { + res: [usize; 2], + bins: usize, +} + +impl CcdBuilder { + pub fn build(&self) -> Array3 { + Array3::zeros([self.res[X], self.res[Y], self.bins]) + } +} \ No newline at end of file diff --git a/src/io/output/mod.rs b/src/io/output/mod.rs index 63ef456..3851b2e 100644 --- a/src/io/output/mod.rs +++ b/src/io/output/mod.rs @@ -1,3 +1,4 @@ +pub mod ccd_builder; pub mod output; pub mod output_config; pub mod output_plane; @@ -10,6 +11,7 @@ pub mod photon_collector; pub mod photon_collector_builder; pub use self::{ + ccd_builder::*, output::*, output_config::*, output_plane::*, diff --git a/src/io/output/output.rs b/src/io/output/output.rs index f87fff1..663d13c 100644 --- a/src/io/output/output.rs +++ b/src/io/output/output.rs @@ -1,7 +1,7 @@ use crate::{ data::Histogram, img::Image, - io::output::{OutputPlane, OutputVolume, OutputRegistry, PhotonCollector}, + io::output::{OutputPlane, OutputRegistry, OutputVolume, PhotonCollector}, phys::Photon, }; use ndarray::Array3; @@ -24,4 +24,21 @@ pub struct Output { /// Contains the mapping between index and name for /// each of the output types. pub reg: OutputRegistry, +} + +impl Output { + /// This function polls each of the output volumes in the output object to + /// find the closest voxel distance based on the position of the current + /// photon packet. This will then return the shortest distance to the + /// the current voxel boundary. There may be a case where there is no voxel + /// in the path of travel of the packet, in that case return `None`. + pub fn voxel_dist(&self, phot: &Photon) -> Option { + let dists = self.vol.iter() + .map(|grid| { grid.voxel_dist(phot) }) + .filter(Option::is_some); + match dists.min_by(|a, b| a.partial_cmp(b).unwrap()) { + Some(val) => val, + None => None, + } + } } \ No newline at end of file diff --git a/src/io/output/output_config.rs b/src/io/output/output_config.rs index 404f87c..d5a95fa 100644 --- a/src/io/output/output_config.rs +++ b/src/io/output/output_config.rs @@ -1,22 +1,118 @@ use std::{ - collections::HashMap, + collections::BTreeMap, fmt }; -use serde::{Serialize}; +use serde::Serialize; use arctk_attr::file; -use crate::{io::output::{OutputPlaneBuilder, OutputVolumeBuilder, PhotonCollectorBuilder}}; +use crate::{ + data::HistogramBuilder, + img::ImageBuilder, + io::output::{OutputPlaneBuilder, OutputVolumeBuilder, PhotonCollectorBuilder, Output}, + ord::Name}; + +use super::{CcdBuilder, OutputRegistry}; #[file] #[derive(Serialize)] pub struct OutputConfig { - pub volumes: Option>, - pub planes: Option>, - pub photon_collectors: Option>, + pub volumes: Option>, + pub planes: Option>, + pub photon_collectors: Option>, + pub spectra: Option>, + pub images: Option>, + pub ccds: Option>, + pub photos: Option>, } impl OutputConfig { + + pub fn build(&self) -> Output { + let reg = OutputRegistry::new_from_config(self); + // Volume output. + let vol = match &self.volumes { + Some(vols) => { + vols.iter().map(|(_key, conf)| { + conf.build() + }) + .collect() + }, + None => vec![] + }; + + let plane = match &self.planes { + Some(planes) => { + planes.iter().map(|(_key, conf)| { + conf.build() + }) + .collect() + }, + None => vec![] + }; + + let phot_cols = match &self.photon_collectors { + Some(pcs) => { + pcs.iter().map(|(_key, conf)| { + conf.build() + }) + .collect() + }, + None => vec![] + }; + + let specs = match &self.spectra { + Some(specs) => { + specs.iter().map(|(_key, conf)| { + conf.build() + }) + .collect() + }, + None => vec![] + }; + + let imgs = match &self.images { + Some(imgs) => { + imgs.iter().map(|(_key, conf)| { + conf.build() + }) + .collect() + }, + None => vec![] + }; + + let ccds = match &self.ccds { + Some(ccds) => { + ccds.iter().map(|(_key, conf)| { + conf.build() + }) + .collect() + }, + None => vec![] + }; + + let photos = match &self.photos { + Some(phots) => { + phots.iter().map(|(_key, conf)| { + conf.build() + }) + .collect() + }, + None => vec![] + }; + + Output { + vol, + plane, + phot_cols, + specs, + imgs, + ccds, + photos, + reg, + } + } + pub fn n_volumes(&self) -> usize { match &self.volumes { Some(vol) => vol.iter().count(), @@ -37,6 +133,34 @@ impl OutputConfig { None => 0, } } + + pub fn n_spectra(&self) -> usize { + match &self.spectra { + Some(spec) => spec.iter().count(), + None => 0, + } + } + + pub fn n_images(&self) -> usize { + match &self.images { + Some(img) => img.iter().count(), + None => 0, + } + } + + pub fn n_ccds(&self) -> usize { + match &self.ccds { + Some(ccd) => ccd.iter().count(), + None => 0, + } + } + + pub fn n_photos(&self) -> usize { + match &self.photos { + Some(phot) => phot.iter().count(), + None => 0, + } + } } impl fmt::Display for OutputConfig { @@ -66,6 +190,11 @@ impl fmt::Display for OutputConfig { } write!(f, " }}\n")?; } + + // TODO: Add output for spectra. + // TODO: Add output for images. + // TODO: Add output for CCDs. + // TODO: Add output for photos. write!(f, "}}") } @@ -91,6 +220,22 @@ mod tests { photon_collectors: { terrain_collector: { kill_photons: false }, sky_collector: { kill_photons: true }, + }, + spectra: { + spec: {min: 400E-9, max: 800E-9, bins: 500}, + }, + images: { + small_image: { res: [1024, 768] }, + larger_image: { res: [1920, 1080] }, + uhd_image: { res: [3840, 2160] }, + }, + ccds: { + default_ccd: { res: [1024, 768], bins: 10}, + }, + photos: { + small_image: { res: [1024, 768] }, + larger_image: { res: [1920, 1080] }, + uhd_image: { res: [3840, 2160] }, } } "#; @@ -102,6 +247,49 @@ mod tests { assert_eq!(conf.n_volumes(), 2); assert_eq!(conf.n_planes(), 1); assert_eq!(conf.n_photon_collectors(), 2); + assert_eq!(conf.n_spectra(), 1); + assert_eq!(conf.n_images(), 3); + assert_eq!(conf.n_ccds(), 1); + assert_eq!(conf.n_photos(), 3); } + #[test] + fn test_build_output() { + let conf_str = r#" + { + volumes: { + full_vol: { boundary: [[0, 0, 0], [10, 10, 10]], res: [8000, 8000, 1000], param: "energy" }, + partial_vol: { boundary: [[2.5, 2.5, 0], [2.5, 2.5, 10]], res: [100, 100, 10], param: "energy" }, + }, + planes: { + bottom: { boundary: [[0, 0], [10, 10]], res: [10, 10], plane: "xy" }, + }, + photon_collectors: { + terrain_collector: { kill_photons: false }, + sky_collector: { kill_photons: true }, + }, + spectra: { + spec: {min: 400E-9, max: 800E-9, bins: 500}, + }, + images: { + small_image: { res: [1024, 768] }, + larger_image: { res: [1920, 1080] }, + uhd_image: { res: [3840, 2160] }, + }, + ccds: { + default_ccd: { res: [1024, 768], bins: 10}, + }, + photos: { + small_image: { res: [1024, 768] }, + larger_image: { res: [1920, 1080] }, + uhd_image: { res: [3840, 2160] }, + } + } + "#; + // Deserialise from the provided string above. + let conf: OutputConfig = json5::from_str(conf_str).unwrap(); + let _out = conf.build(); + + // TODO: Implement some tests for the output object here. + } } \ No newline at end of file diff --git a/src/io/output/output_registry.rs b/src/io/output/output_registry.rs index c16efd9..9504944 100644 --- a/src/io/output/output_registry.rs +++ b/src/io/output/output_registry.rs @@ -1,6 +1,5 @@ use crate::{ - ord::Register, - io::output::Output, + io::output::OutputConfig, ord::{Name, Register} }; pub struct OutputRegistry { @@ -13,5 +12,39 @@ pub struct OutputRegistry { } impl OutputRegistry { - + pub fn new_from_config(out: &OutputConfig) -> OutputRegistry { + // Get the keys for the volume outputs. + let vol_keys = match &out.volumes { + Some(vols) => vols.keys().map(|k| k.clone()).collect(), + None => vec![] + }; + let vol_reg = Register::new(vol_keys); + + // Get the keys for the plane outputs. + let plane_keys = match &out.planes { + Some(planes) => planes.keys().map(|k| k.clone()).collect(), + None => vec![] + }; + let plane_reg = Register::new(plane_keys); + + // Get the keys for photon collector outputs. + let photcol_keys = match &out.photon_collectors { + Some(photcol) => photcol.keys().map(|k| k.clone()).collect(), + None => vec![] + }; + let phot_cols_reg = Register::new(photcol_keys); + + let spec_reg = Register::new(vec![]); + let img_reg = Register::new(vec![]); + let ccd_reg = Register::new(vec![]); + + Self { + vol_reg, + plane_reg, + phot_cols_reg, + spec_reg, + img_reg, + ccd_reg, + } + } } \ No newline at end of file diff --git a/src/io/output/output_volume.rs b/src/io/output/output_volume.rs index c81ed72..7b648e7 100644 --- a/src/io/output/output_volume.rs +++ b/src/io/output/output_volume.rs @@ -5,7 +5,12 @@ use pyo3::prelude::*; use ndarray::Array3; use serde::{Deserialize, Serialize}; use crate::{ - access, core::Real, fs::Save, geom::Cube, ord::cartesian::{X, Y, Z} + access, + fs::Save, + geom::{Cube, Trace}, + math::{Point3, Vec3}, + ord::cartesian::{X, Y, Z}, + phys::Photon }; #[derive(PartialEq, Debug, Clone, Serialize, Deserialize)] @@ -22,6 +27,7 @@ pub struct OutputVolume { boundary: Cube, res: [usize; 3], param: OutputParameter, + voxel_size: Vec3, data: Array3, } @@ -36,17 +42,74 @@ impl OutputVolume { debug_assert!(res[Y] > 0); debug_assert!(res[Z] > 0); + let mut voxel_size = boundary.widths(); + for (w, n) in voxel_size.iter_mut().zip(&res) { + *w /= *n as f64; + } + Self { boundary, res, param, + voxel_size, data: Array3::zeros(res) } } + #[inline] + #[must_use] pub fn cell_volume(&self) -> f64 { self.boundary.vol() / (self.res[X] * self.res[Y] * self.res[Z]) as f64 } + + /// If the given position is contained within the grid, + /// generate the index for the given position within the grid. + #[inline] + #[must_use] + pub fn gen_index(&self, p: &Point3) -> Option<[usize; 3]> { + self.boundary.contains(p).then(|| { + let mins = self.boundary.mins(); + let maxs = self.boundary.maxs(); + + [ + (((p.x() - mins.x()) / (maxs.x() - mins.x())) * self.res[X] as f64).floor() + as usize, + (((p.y() - mins.y()) / (maxs.y() - mins.y())) * self.res[Y] as f64).floor() + as usize, + (((p.z() - mins.z()) / (maxs.z() - mins.z())) * self.res[Z] as f64).floor() + as usize, + ] + }) + } + + /// If the given position is contained within the grid, + /// generate the index and voxel for the given position within the grid. + #[inline] + #[must_use] + pub fn gen_index_voxel(&self, p: &Point3) -> Option<([usize; 3], Cube)> { + if let Some(index) = self.gen_index(p) { + let mut min = *self.boundary.mins(); + *min.x_mut() += self.voxel_size[X] * index[X] as f64; + *min.y_mut() += self.voxel_size[Y] * index[Y] as f64; + *min.z_mut() += self.voxel_size[Z] * index[Z] as f64; + + let boundary = Cube::new(min, min + self.voxel_size); + debug_assert!(boundary.contains(p)); + + Some((index, boundary)) + } else { + None + } + } + + /// Returns the distance to the nearest voxel boundary, if one exists. + #[inline] + #[must_use] + pub fn voxel_dist(&self, phot: &Photon) -> Option { + let (_index, voxel) = self.gen_index_voxel(phot.ray().pos())?; + let dist = voxel.dist(phot.ray())?; + Some(dist) + } } impl AddAssign<&Self> for OutputVolume { @@ -70,7 +133,6 @@ impl Save for OutputVolume { #[cfg(test)] mod tests { use crate::math::Point3; - use json5; use super::*; #[test] From f57a8bec2ad4374122d3fa737454d377a281beff Mon Sep 17 00:00:00 2001 From: Sam Morrell Date: Wed, 14 Aug 2024 01:30:37 +0100 Subject: [PATCH 09/16] Added `output` section to parameters input. This now means that the output parameters can be read in as part of the main JSON5 input for the MCRT software. I had added a redirect here, so it can be included either in the same file, or separate. Also added output section to display. --- src/sim/param/parameters.rs | 11 ++++++----- src/sim/param/parameters_builder.rs | 14 ++++++++------ src/sim/param/parameters_builder_loader.rs | 14 ++++++-------- 3 files changed, 20 insertions(+), 19 deletions(-) diff --git a/src/sim/param/parameters.rs b/src/sim/param/parameters.rs index 352c000..07bda3a 100644 --- a/src/sim/param/parameters.rs +++ b/src/sim/param/parameters.rs @@ -1,11 +1,7 @@ //! Runtime parameters. use crate::{ - fmt_report, - geom::{Grid, SurfaceLinker, TreeSettings}, - ord::Set, - phys::{LightLinker, Material}, - sim::{AttributeLinkerLinkerLinkerLinkerLinker, Engine, Settings}, + fmt_report, geom::{Grid, SurfaceLinker, TreeSettings}, io::output::OutputConfig, ord::Set, phys::{LightLinker, Material}, sim::{AttributeLinkerLinkerLinkerLinkerLinker, Engine, Settings} }; use std::fmt::{Display, Error, Formatter}; @@ -27,6 +23,8 @@ pub struct Parameters { pub lights: Set, /// Engine selection. pub engine: Engine, + /// Outputs + pub output: OutputConfig, } impl Parameters { @@ -43,6 +41,7 @@ impl Parameters { mats: Set, lights: Set, engine: Engine, + output: OutputConfig, ) -> Self { Self { sett, @@ -53,6 +52,7 @@ impl Parameters { mats, lights, engine, + output, } } } @@ -69,6 +69,7 @@ impl Display for Parameters { fmt_report!(fmt, self.mats, "materials"); fmt_report!(fmt, self.lights, "lights"); fmt_report!(fmt, self.engine, "engine"); + fmt_report!(fmt, self.output, "output"); Ok(()) } } diff --git a/src/sim/param/parameters_builder.rs b/src/sim/param/parameters_builder.rs index d8cf669..449fa02 100644 --- a/src/sim/param/parameters_builder.rs +++ b/src/sim/param/parameters_builder.rs @@ -1,11 +1,7 @@ //! Buildable parameters. use crate::{ - fmt_report, - geom::{GridBuilder, SurfaceLinker, TreeSettings}, - ord::{Build, Set}, - phys::{LightLinkerBuilder, MaterialBuilder}, - sim::{AttributeLinkerLinkerLinkerLinkerLinker, EngineBuilder, Parameters, Settings}, + fmt_report, geom::{GridBuilder, SurfaceLinker, TreeSettings}, io::output::OutputConfig, ord::{Build, Set}, phys::{LightLinkerBuilder, MaterialBuilder}, sim::{AttributeLinkerLinkerLinkerLinkerLinker, EngineBuilder, Parameters, Settings} }; use std::fmt::{Display, Error, Formatter}; @@ -27,6 +23,8 @@ pub struct ParametersBuilder { lights: Set, /// Engine selection. engine: EngineBuilder, + /// Output + output: OutputConfig, } impl ParametersBuilder { @@ -43,6 +41,7 @@ impl ParametersBuilder { mats: Set, lights: Set, engine: EngineBuilder, + output: OutputConfig ) -> Self { Self { sett, @@ -53,6 +52,7 @@ impl ParametersBuilder { mats, lights, engine, + output, } } } @@ -70,8 +70,9 @@ impl Build for ParametersBuilder { let mats = self.mats.build(); let light = self.lights.build(); let engine = self.engine.build(); + let output = self.output; - Self::Inst::new(sett, tree, grid, surfs, attrs, mats, light, engine) + Self::Inst::new(sett, tree, grid, surfs, attrs, mats, light, engine, output) } } @@ -87,6 +88,7 @@ impl Display for ParametersBuilder { fmt_report!(fmt, self.mats, "materials"); fmt_report!(fmt, self.lights, "lights"); fmt_report!(fmt, self.engine, "engine"); + fmt_report!(fmt, self.output, "output"); Ok(()) } } diff --git a/src/sim/param/parameters_builder_loader.rs b/src/sim/param/parameters_builder_loader.rs index f1afaea..596d1e6 100644 --- a/src/sim/param/parameters_builder_loader.rs +++ b/src/sim/param/parameters_builder_loader.rs @@ -1,14 +1,9 @@ //! Loadable parameters. use crate::{ - err::Error, - fs::{Load, Redirect}, - geom::{GridBuilder, SurfaceLinkerLoader, TreeSettings}, - ord::Set, - phys::{LightLinkerBuilderLoader, MaterialBuilder}, - sim::{ + err::Error, fs::{Load, Redirect}, geom::{GridBuilder, SurfaceLinkerLoader, TreeSettings}, io::output::OutputConfig, ord::Set, phys::{LightLinkerBuilderLoader, MaterialBuilder}, sim::{ AttributeLinkerLinkerLinkerLinkerLinker, EngineBuilderLoader, ParametersBuilder, Settings, - }, + } }; use arctk_attr::file; use std::path::Path; @@ -32,6 +27,8 @@ pub struct ParametersBuilderLoader { lights: Redirect>, /// Engine selection. engine: EngineBuilderLoader, + /// Output + output: Redirect, } impl Load for ParametersBuilderLoader { @@ -47,9 +44,10 @@ impl Load for ParametersBuilderLoader { let mats = self.mats.load(in_dir)?.load(in_dir)?; let lights = self.lights.load(in_dir)?.load(in_dir)?; let engine = self.engine.load(in_dir)?; + let output = self.output.load(in_dir)?; Ok(Self::Inst::new( - sett, tree, grid, surfs, attrs, mats, lights, engine, + sett, tree, grid, surfs, attrs, mats, lights, engine, output, )) } } From c59df61eefc6b80395bc3f9d0e831a765b3ded63 Mon Sep 17 00:00:00 2001 From: Sam Morrell Date: Wed, 14 Aug 2024 01:54:19 +0100 Subject: [PATCH 10/16] Fixed `Output` voxel_dist function. This function was previously using an invalid min finding iter-chain. It also now sensibly returns f64::INFINITY instead of a None in the case of distance to a voxel being None (no voxel boundaries in path of packet). --- src/io/output/output.rs | 12 +++++++----- 1 file changed, 7 insertions(+), 5 deletions(-) diff --git a/src/io/output/output.rs b/src/io/output/output.rs index 663d13c..4c64d6e 100644 --- a/src/io/output/output.rs +++ b/src/io/output/output.rs @@ -32,13 +32,15 @@ impl Output { /// photon packet. This will then return the shortest distance to the /// the current voxel boundary. There may be a case where there is no voxel /// in the path of travel of the packet, in that case return `None`. - pub fn voxel_dist(&self, phot: &Photon) -> Option { - let dists = self.vol.iter() + pub fn voxel_dist(&self, phot: &Photon) -> f64 { + let dists: Vec = self.vol.iter() .map(|grid| { grid.voxel_dist(phot) }) - .filter(Option::is_some); - match dists.min_by(|a, b| a.partial_cmp(b).unwrap()) { + .filter(Option::is_some) + .map(Option::unwrap) + .collect(); + match dists.into_iter().reduce(f64::min) { Some(val) => val, - None => None, + None => f64::INFINITY, } } } \ No newline at end of file From 64e8a8aa557cecf1b781651e31038d9c6e394c71 Mon Sep 17 00:00:00 2001 From: Sam Morrell Date: Wed, 14 Aug 2024 12:53:54 +0100 Subject: [PATCH 11/16] Starting integration of new `Output` object into MCRT This commit adds new the new Output struct to the MCRT executable. Including implementing new types throughout the code and integrating into existing code. Some more specific comments on changes: - Modified engines and `travel` implementation to output into new objects. - Modified `mcrt` exec to use new output object. - Added extra implementations to clone output. - Added placeholder implementations of `Save` for `Output`. --- src/bin/mcrt.rs | 214 ++++++++++++------------- src/geom/domain/boundary.rs | 6 +- src/geom/domain/grid.rs | 3 - src/io/output/output.rs | 64 +++++++- src/io/output/output_plane.rs | 12 +- src/io/output/output_registry.rs | 1 + src/io/output/output_volume.rs | 51 +++++- src/io/output/output_volume_builder.rs | 2 +- src/ord/register.rs | 1 + src/ord/set.rs | 2 +- src/sim/engine/engine.rs | 3 +- src/sim/engine/engines/fluorescence.rs | 23 +-- src/sim/engine/engines/photo.rs | 23 +-- src/sim/engine/engines/raman.rs | 23 +-- src/sim/engine/engines/standard.rs | 33 ++-- src/sim/run.rs | 13 +- src/sim/surface.rs | 3 +- src/sim/travel.rs | 33 ++-- 18 files changed, 324 insertions(+), 186 deletions(-) diff --git a/src/bin/mcrt.rs b/src/bin/mcrt.rs index 39feae2..88a38df 100644 --- a/src/bin/mcrt.rs +++ b/src/bin/mcrt.rs @@ -1,22 +1,19 @@ //! Monte-Carlo radiative transfer simulation binary. //! Compute the radiative field for a given set of setup and light source. -use ndarray::Array3; use std::{ env::current_dir, path::{Path, PathBuf}, }; use aetherus::{ args, - data::Histogram, fs::{File, Load, Save}, - geom::{Grid, Tree, Boundary}, - img::{Colour, Image}, - ord::{Build, Link, Register, Set, cartesian::{X, Y}}, + geom::{Tree, Boundary}, + ord::{Build, Link, Register, Set}, report, - io::output::PhotonCollector, + io::output::{Output, PhotonCollector}, sim::{ - run, AttributeLinkerLinkerLinkerLinkerLinker as Attr, Engine, Input, Output, Parameters, + run, AttributeLinkerLinkerLinkerLinkerLinker as Attr, Engine, Input, Parameters, ParametersBuilderLoader, }, util::{ @@ -25,6 +22,7 @@ use aetherus::{ fmt::term, }, }; +use nalgebra::base; /// Backup print width if the terminal width can not be determined. const BACKUP_TERM_WIDTH: usize = 80; @@ -50,15 +48,17 @@ fn main() { sub_section(term_width, "Registration"); let (spec_reg, img_reg, ccd_reg, phot_col_reg) = gen_detector_registers(¶ms.attrs); - let base_output = gen_base_output( - &engine, - &grid, - &spec_reg, - &img_reg, - &ccd_reg, - &phot_col_reg, - ¶ms.attrs, - ); + // let base_output = gen_base_output( + // &engine, + // &grid, + // &spec_reg, + // &img_reg, + // &ccd_reg, + // &phot_col_reg, + // ¶ms.attrs, + // ); + + let base_output = params.output.build(); sub_section(term_width, "Linking"); let lights = params @@ -68,13 +68,13 @@ fn main() { report!(lights, "lights"); let attrs = params .attrs - .link(phot_col_reg.set()) + .link(base_output.reg.phot_cols_reg.set()) .expect("Failed to link photon collectors to attributes.") - .link(ccd_reg.set()) + .link(base_output.reg.ccd_reg.set()) .expect("Failed to link ccds to attributes.") - .link(img_reg.set()) + .link(base_output.reg.img_reg.set()) .expect("Failed to link imagers to attributes.") - .link(spec_reg.set()) + .link(base_output.reg.spec_reg.set()) .expect("Failed to link spectrometers to attributes.") .link(&mats) .expect("Failed to link materials to attributes."); @@ -127,12 +127,12 @@ fn main() { } } - output += &data; + output += data; output }); section(term_width, "Saving"); - report!(data, "data"); + //report!(data, "data"); data.save(&out_dir).expect("Failed to save output data."); section(term_width, "Finished"); @@ -212,88 +212,88 @@ fn gen_detector_registers(attrs: &Set) -> (Register, Register, Register, R (spec_reg, img_reg, ccd_reg, phot_col_reg) } -/// Generate the base output instance. -fn gen_base_output<'a>( - engine: &Engine, - grid: &Grid, - spec_reg: &'a Register, - img_reg: &'a Register, - ccd_reg: &'a Register, - phot_col_reg: &'a Register, - attrs: &Set, -) -> Output<'a> { - let res = *grid.res(); - - let mut specs = Vec::with_capacity(spec_reg.len()); - for name in spec_reg.set().map().keys() { - for attr in attrs.values() { - if let Attr::Spectrometer(spec_name, [min, max], bins) = attr { - if name == spec_name { - specs.push(Histogram::new(*min, *max, *bins)); - continue; - } - } - } - } - - let mut imgs = Vec::with_capacity(img_reg.len()); - let background = Colour::new(0.0, 0.0, 0.0, 1.0); - for name in img_reg.set().map().keys() { - for attr in attrs.values() { - if let Attr::Imager(img_name, res, _width, _center, _forward) = attr { - if name == img_name { - imgs.push(Image::new_blank(*res, background)); - continue; - } - } - } - } - - let mut ccds = Vec::with_capacity(ccd_reg.len()); - for name in ccd_reg.set().map().keys() { - for attr in attrs.values() { - if let Attr::Ccd(ccd_name, res, _width, _center, _forward, binner) = attr { - if name == ccd_name { - ccds.push(Array3::zeros([res[X], res[Y], binner.bins() as usize])); - continue; - } - } - } - } - - let mut photos = Vec::new(); - if let Engine::Photo(frames, res) = engine { - photos.reserve(frames.len()); - for _ in 0..frames.len() { - photos.push(Image::new_blank(*res, background)); - } - } - - let mut phot_cols: Vec = Vec::new(); - for name in phot_col_reg.set().map().keys() { - for attr in attrs.values() { - if let Attr::PhotonCollector(phot_col_id, kill_phot) = attr { - if name == phot_col_id { - let mut photcol = PhotonCollector::new(); - photcol.kill_photon = *kill_phot; - phot_cols.push(photcol); - continue; - } - } - } - } - - Output::new( - grid.boundary().clone(), - res, - spec_reg, - img_reg, - ccd_reg, - phot_col_reg, - specs, - imgs, - ccds, - photos, - phot_cols, - ) -} +// Generate the base output instance. +// fn gen_base_output<'a>( +// engine: &Engine, +// grid: &Grid, +// spec_reg: &'a Register, +// img_reg: &'a Register, +// ccd_reg: &'a Register, +// phot_col_reg: &'a Register, +// attrs: &Set, +// ) -> Output<'a> { +// let res = *grid.res(); + +// let mut specs = Vec::with_capacity(spec_reg.len()); +// for name in spec_reg.set().map().keys() { +// for attr in attrs.values() { +// if let Attr::Spectrometer(spec_name, [min, max], bins) = attr { +// if name == spec_name { +// specs.push(Histogram::new(*min, *max, *bins)); +// continue; +// } +// } +// } +// } + +// let mut imgs = Vec::with_capacity(img_reg.len()); +// let background = Colour::new(0.0, 0.0, 0.0, 1.0); +// for name in img_reg.set().map().keys() { +// for attr in attrs.values() { +// if let Attr::Imager(img_name, res, _width, _center, _forward) = attr { +// if name == img_name { +// imgs.push(Image::new_blank(*res, background)); +// continue; +// } +// } +// } +// } + +// let mut ccds = Vec::with_capacity(ccd_reg.len()); +// for name in ccd_reg.set().map().keys() { +// for attr in attrs.values() { +// if let Attr::Ccd(ccd_name, res, _width, _center, _forward, binner) = attr { +// if name == ccd_name { +// ccds.push(Array3::zeros([res[X], res[Y], binner.bins() as usize])); +// continue; +// } +// } +// } +// } + +// let mut photos = Vec::new(); +// if let Engine::Photo(frames, res) = engine { +// photos.reserve(frames.len()); +// for _ in 0..frames.len() { +// photos.push(Image::new_blank(*res, background)); +// } +// } + +// let mut phot_cols: Vec = Vec::new(); +// for name in phot_col_reg.set().map().keys() { +// for attr in attrs.values() { +// if let Attr::PhotonCollector(phot_col_id, kill_phot) = attr { +// if name == phot_col_id { +// let mut photcol = PhotonCollector::new(); +// photcol.kill_photon = *kill_phot; +// phot_cols.push(photcol); +// continue; +// } +// } +// } +// } + +// Output::new( +// grid.boundary().clone(), +// res, +// spec_reg, +// img_reg, +// ccd_reg, +// phot_col_reg, +// specs, +// imgs, +// ccds, +// photos, +// phot_cols, +// ) +// } diff --git a/src/geom/domain/boundary.rs b/src/geom/domain/boundary.rs index 8eaca6c..ddad49a 100644 --- a/src/geom/domain/boundary.rs +++ b/src/geom/domain/boundary.rs @@ -1,5 +1,9 @@ use crate::{ - access, clone, fmt_report, geom::{plane::ray_plane_intersection, Cube, Hit, Ray, Side, Trace}, math::{Dir3, Point3, Vec3}, phys::{Photon, Reflectance}, sim::Attribute + access, clone, fmt_report, + geom::{plane::ray_plane_intersection, Cube, Hit, Ray, Side, Trace}, + math::{Dir3, Point3, Vec3}, + phys::{Photon, Reflectance}, + sim::Attribute }; use rand::rngs::ThreadRng; use std::fmt::{Display, Formatter}; diff --git a/src/geom/domain/grid.rs b/src/geom/domain/grid.rs index 04e7712..f57702b 100644 --- a/src/geom/domain/grid.rs +++ b/src/geom/domain/grid.rs @@ -152,9 +152,6 @@ impl Display for Grid { #[cfg(test)] mod tests { use super::*; - use crate::{ - geom::Cube, - }; #[test] fn test_new() { diff --git a/src/io/output/output.rs b/src/io/output/output.rs index 4c64d6e..1c60226 100644 --- a/src/io/output/output.rs +++ b/src/io/output/output.rs @@ -1,10 +1,19 @@ use crate::{ + fs::Save, data::Histogram, img::Image, + err::Error, io::output::{OutputPlane, OutputRegistry, OutputVolume, PhotonCollector}, phys::Photon, }; +use std::{ + ops::AddAssign, + path::Path, +}; use ndarray::Array3; +use super::OutputParameter; + +#[derive(Clone)] pub struct Output { /// Output volumes. pub vol: Vec, @@ -43,4 +52,57 @@ impl Output { None => f64::INFINITY, } } -} \ No newline at end of file + + pub fn get_volumes_for_param(&self, param: OutputParameter) -> Vec<&OutputVolume> { + self.vol.iter() + .filter(|&vol| vol.param() == ¶m) + .collect() + } + + pub fn get_volumes_for_param_mut(&mut self, param: OutputParameter) -> Vec<&mut OutputVolume> { + self.vol.iter_mut() + .filter(|vol| vol.param() == ¶m) + .collect() + } +} + +impl AddAssign for Output { + #[inline] + fn add_assign(&mut self, rhs: Self) { + for (a, b) in self.vol.iter_mut().zip(&rhs.vol) { + *a += b; + } + + for (a, b) in self.plane.iter_mut().zip(&rhs.plane) { + *a += b; + } + + for (a, b) in self.phot_cols.iter_mut().zip(&rhs.phot_cols) { + *a += b; + } + + for (a, b) in self.specs.iter_mut().zip(&rhs.specs) { + *a += b; + } + + for (a, b) in self.imgs.iter_mut().zip(&rhs.imgs) { + *a += b; + } + + for (a, b) in self.ccds.iter_mut().zip(&rhs.ccds) { + *a += b; + } + + for (a, b) in self.photos.iter_mut().zip(&rhs.photos) { + *a += b; + } + } +} + + +impl Save for Output { + #[inline] + fn save_data(&self, out_dir: &Path) -> Result<(), Error> { + todo!() + } +} diff --git a/src/io/output/output_plane.rs b/src/io/output/output_plane.rs index 06b59da..8d16fc2 100644 --- a/src/io/output/output_plane.rs +++ b/src/io/output/output_plane.rs @@ -1,3 +1,5 @@ +use std::ops::AddAssign; + use ndarray::Array2; use serde::{Deserialize, Serialize}; use crate::{ @@ -24,7 +26,7 @@ impl AxisAlignedPlane { } } -#[derive(Debug)] +#[derive(Debug, Clone)] pub struct OutputPlane { mins: Point2, maxs: Point2, @@ -96,6 +98,14 @@ impl OutputPlane { } +impl AddAssign<&Self> for OutputPlane { + fn add_assign(&mut self, rhs: &Self) { + debug_assert_eq!(self.res, rhs.res); + + self.data += &rhs.data; + } +} + #[cfg(test)] mod tests { use super::*; diff --git a/src/io/output/output_registry.rs b/src/io/output/output_registry.rs index 9504944..6e0a57b 100644 --- a/src/io/output/output_registry.rs +++ b/src/io/output/output_registry.rs @@ -2,6 +2,7 @@ use crate::{ io::output::OutputConfig, ord::{Name, Register} }; +#[derive(Clone)] pub struct OutputRegistry { pub vol_reg: Register, pub plane_reg: Register, diff --git a/src/io/output/output_volume.rs b/src/io/output/output_volume.rs index 7b648e7..ca5fed5 100644 --- a/src/io/output/output_volume.rs +++ b/src/io/output/output_volume.rs @@ -2,10 +2,12 @@ use std::ops::AddAssign; #[cfg(feature = "pyo3")] use pyo3::prelude::*; +use std::fmt::{Display, Formatter}; use ndarray::Array3; use serde::{Deserialize, Serialize}; use crate::{ access, + fmt_report, fs::Save, geom::{Cube, Trace}, math::{Point3, Vec3}, @@ -17,12 +19,13 @@ use crate::{ #[cfg_attr(feature = "pyo3", pyclass)] #[serde(rename_all = "lowercase")] pub enum OutputParameter { + Emission, Energy, Absorption, Shift, } -#[derive(Debug)] +#[derive(Debug, Clone)] pub struct OutputVolume { boundary: Cube, res: [usize; 3], @@ -34,6 +37,7 @@ pub struct OutputVolume { impl OutputVolume { access!(boundary: Cube); access!(res: [usize; 3]); + access!(param: OutputParameter); access!(data, data_mut: Array3); pub fn new(boundary: Cube, res: [usize; 3], param: OutputParameter) -> Self { @@ -58,10 +62,23 @@ impl OutputVolume { #[inline] #[must_use] - pub fn cell_volume(&self) -> f64 { + pub fn voxel_volume(&self) -> f64 { self.boundary.vol() / (self.res[X] * self.res[Y] * self.res[Z]) as f64 } + /// Determine the total number of cells. + #[inline] + #[must_use] + pub const fn num_cells(&self) -> usize { + self.res[X] * self.res[Y] * self.res[Z] + } + + #[inline] + #[must_use] + pub fn contains(&self, p: &Point3) -> bool { + self.boundary.contains(p) + } + /// If the given position is contained within the grid, /// generate the index for the given position within the grid. #[inline] @@ -115,7 +132,7 @@ impl OutputVolume { impl AddAssign<&Self> for OutputVolume { fn add_assign(&mut self, rhs: &Self) { debug_assert_eq!(self.boundary(), rhs.boundary()); - debug_assert_eq!(self.cell_volume(), rhs.cell_volume()); + debug_assert_eq!(self.voxel_volume(), rhs.voxel_volume()); debug_assert_eq!(self.param, rhs.param); self.data += &rhs.data @@ -125,7 +142,31 @@ impl AddAssign<&Self> for OutputVolume { impl Save for OutputVolume { #[inline] fn save_data(&self, path: &std::path::Path) -> Result<(), crate::err::Error> { - (&self.data / self.cell_volume()).save(&path)?; + (&self.data / self.voxel_volume()).save(&path)?; + Ok(()) + } +} + +impl Display for OutputVolume { + #[inline] + fn fmt(&self, fmt: &mut Formatter) -> Result<(), std::fmt::Error> { + writeln!(fmt, "...")?; + fmt_report!(fmt, self.boundary, "boundary"); + fmt_report!( + fmt, + &format!("[{} x {} x {}]", self.res[X], self.res[Y], self.res[Z]), + "resolution" + ); + fmt_report!( + fmt, + &format!( + "({}, {}, {})", + self.voxel_size.x(), + self.voxel_size.y(), + self.voxel_size.z() + ), + "voxel size" + ); Ok(()) } } @@ -139,6 +180,6 @@ mod tests { fn new_test() { let boundary = Cube::new(Point3::new(0.0,0.0,0.0), Point3::new(10.0,10.0,10.0)); let outvol = OutputVolume::new(boundary, [10, 10, 10], OutputParameter::Energy); - assert_eq!(outvol.cell_volume(), 1.0); + assert_eq!(outvol.voxel_volume(), 1.0); } } diff --git a/src/io/output/output_volume_builder.rs b/src/io/output/output_volume_builder.rs index 184b3fc..af6b234 100644 --- a/src/io/output/output_volume_builder.rs +++ b/src/io/output/output_volume_builder.rs @@ -38,6 +38,6 @@ mod tests { let builder: OutputVolumeBuilder = json5::from_str(input).unwrap(); let outvol = builder.build(); - assert_eq!(outvol.cell_volume(), 1.0); + assert_eq!(outvol.voxel_volume(), 1.0); } } diff --git a/src/ord/register.rs b/src/ord/register.rs index 5b60725..941b3fe 100644 --- a/src/ord/register.rs +++ b/src/ord/register.rs @@ -7,6 +7,7 @@ use crate::{ use std::fmt::{Display, Error, Formatter}; /// Register used to index named data. +#[derive(Clone)] pub struct Register(Set); impl Register { diff --git a/src/ord/set.rs b/src/ord/set.rs index a7fb162..494e98f 100644 --- a/src/ord/set.rs +++ b/src/ord/set.rs @@ -14,7 +14,7 @@ use std::{ }; /// Data map. -#[derive(Debug, Deserialize, Serialize)] +#[derive(Debug, Deserialize, Serialize, Clone)] pub struct Set(Map); impl Set { diff --git a/src/sim/engine/engine.rs b/src/sim/engine/engine.rs index f5f2fd8..158e218 100644 --- a/src/sim/engine/engine.rs +++ b/src/sim/engine/engine.rs @@ -4,7 +4,8 @@ use crate::{ math::{Formula, Point3}, ord::cartesian::{X, Y}, phys::Photon, - sim::{engines, Frame, Input, Output}, + sim::{engines, Frame, Input}, + io::output::Output, }; use ndarray::Array3; use rand::rngs::ThreadRng; diff --git a/src/sim/engine/engines/fluorescence.rs b/src/sim/engine/engines/fluorescence.rs index b571193..a976aa1 100644 --- a/src/sim/engine/engines/fluorescence.rs +++ b/src/sim/engine/engines/fluorescence.rs @@ -4,7 +4,8 @@ use crate::{ geom::Trace, math::Formula, phys::{Local, Photon}, - sim::{scatter::scatter, surface::surface, travel::travel, Event, Input, Output}, + sim::{scatter::scatter, surface::surface, travel::travel, Event, Input}, + io::output::{Output, OutputParameter}, }; use ndarray::Array3; use rand::{rngs::ThreadRng, Rng}; @@ -20,11 +21,11 @@ pub fn fluorescence( mut rng: &mut ThreadRng, mut phot: Photon, ) { - // Check photon is within the grid. - if let Some(index) = input.grid.gen_index(phot.ray().pos()) { - data.emission[index] += phot.power() * phot.weight(); - } else { - panic!("Photon was not emitted within the grid."); + // Add to the emission variables in which the photon is present. + for vol in data.get_volumes_for_param_mut(OutputParameter::Emission) { + if let Some(index) = vol.gen_index(phot.ray().pos()) { + vol.data_mut()[index] += phot.power() * phot.weight(); + } } // Common constants. @@ -84,18 +85,18 @@ pub fn fluorescence( // Event handling. match Event::new(voxel_dist, scat_dist, surf_hit, boundary_hit, bump_dist) { - Event::Voxel(dist) => travel(&mut data, &mut phot, &env, &index, dist + bump_dist), + Event::Voxel(dist) => travel(&mut data, &mut phot, &env, dist + bump_dist), Event::Scattering(dist) => { - travel(&mut data, &mut phot, &env, &index, dist); + travel(&mut data, &mut phot, &env, dist); scatter(&mut rng, &mut phot, &env); } Event::Surface(hit) => { - travel(&mut data, &mut phot, &env, &index, hit.dist()); + travel(&mut data, &mut phot, &env, hit.dist()); surface(&mut rng, &hit, &mut phot, &mut local, &mut data); - travel(&mut data, &mut phot, &env, &index, bump_dist); + travel(&mut data, &mut phot, &env, bump_dist); }, Event::Boundary(boundary_hit) => { - travel(&mut data, &mut phot, &env, &index, boundary_hit.dist()); + travel(&mut data, &mut phot, &env, boundary_hit.dist()); input.bound.apply(rng, &boundary_hit, &mut phot); } } diff --git a/src/sim/engine/engines/photo.rs b/src/sim/engine/engines/photo.rs index 2cc3014..2eb9368 100644 --- a/src/sim/engine/engines/photo.rs +++ b/src/sim/engine/engines/photo.rs @@ -4,9 +4,10 @@ use crate::{ geom::Trace, img::Colour, phys::Photon, + io::output::{Output, OutputParameter}, sim::{ peel_off::peel_off, scatter::scatter, surface::surface, travel::travel, Event, Frame, - Input, Output, + Input, }, }; use rand::{rngs::ThreadRng, Rng}; @@ -21,11 +22,11 @@ pub fn photo( mut rng: &mut ThreadRng, mut phot: Photon, ) { - // Check photon is within the grid. - if let Some(index) = input.grid.gen_index(phot.ray().pos()) { - data.emission[index] += phot.power() * phot.weight(); - } else { - panic!("Photon was not emitted within the grid."); + // Add to the emission variables in which the photon is present. + for vol in data.get_volumes_for_param_mut(OutputParameter::Emission) { + if let Some(index) = vol.gen_index(phot.ray().pos()) { + vol.data_mut()[index] += phot.power() * phot.weight(); + } } // Common constants. @@ -75,9 +76,9 @@ pub fn photo( // Event handling. match Event::new(voxel_dist, scat_dist, surf_hit, boundary_hit, bump_dist) { - Event::Voxel(dist) => travel(&mut data, &mut phot, &env, &index, dist + bump_dist), + Event::Voxel(dist) => travel(&mut data, &mut phot, &env, dist + bump_dist), Event::Scattering(dist) => { - travel(&mut data, &mut phot, &env, &index, dist); + travel(&mut data, &mut phot, &env, dist); // Capture. for (frame, photo) in frames.iter().zip(data.photos.iter_mut()) { @@ -97,12 +98,12 @@ pub fn photo( scatter(&mut rng, &mut phot, &env); } Event::Surface(hit) => { - travel(&mut data, &mut phot, &env, &index, hit.dist()); + travel(&mut data, &mut phot, &env, hit.dist()); surface(&mut rng, &hit, &mut phot, &mut env, &mut data); - travel(&mut data, &mut phot, &env, &index, bump_dist); + travel(&mut data, &mut phot, &env, bump_dist); }, Event::Boundary(boundary_hit) => { - travel(&mut data, &mut phot, &env, &index, boundary_hit.dist()); + travel(&mut data, &mut phot, &env, boundary_hit.dist()); input.bound.apply(rng, &boundary_hit, &mut phot); } } diff --git a/src/sim/engine/engines/raman.rs b/src/sim/engine/engines/raman.rs index de1431f..48244fd 100644 --- a/src/sim/engine/engines/raman.rs +++ b/src/sim/engine/engines/raman.rs @@ -4,7 +4,8 @@ use crate::{ geom::Trace, math::Point3, phys::Photon, - sim::{scatter::shift_scatter, surface::surface, travel::travel, Event, Input, Output}, + sim::{scatter::shift_scatter, surface::surface, travel::travel, Event, Input}, + io::output::{Output, OutputParameter}, }; use rand::{rngs::ThreadRng, Rng}; @@ -18,11 +19,11 @@ pub fn raman( mut rng: &mut ThreadRng, mut phot: Photon, ) { - // Check photon is within the grid. - if let Some(index) = input.grid.gen_index(phot.ray().pos()) { - data.emission[index] += phot.power() * phot.weight(); - } else { - panic!("Photon was not emitted within the grid."); + // Add to the emission variables in which the photon is present. + for vol in data.get_volumes_for_param_mut(OutputParameter::Emission) { + if let Some(index) = vol.gen_index(phot.ray().pos()) { + vol.data_mut()[index] += phot.power() * phot.weight(); + } } // Common constants. @@ -72,9 +73,9 @@ pub fn raman( // Event handling. match Event::new(voxel_dist, scat_dist, surf_hit, boundary_hit, bump_dist) { - Event::Voxel(dist) => travel(&mut data, &mut phot, &env, &index, dist + bump_dist), + Event::Voxel(dist) => travel(&mut data, &mut phot, &env, dist + bump_dist), Event::Scattering(dist) => { - travel(&mut data, &mut phot, &env, &index, dist); + travel(&mut data, &mut phot, &env, dist); // // Capture. // if let Some(weight) = @@ -86,12 +87,12 @@ pub fn raman( shift_scatter(&mut rng, &mut phot, &env); } Event::Surface(hit) => { - travel(&mut data, &mut phot, &env, &index, hit.dist()); + travel(&mut data, &mut phot, &env, hit.dist()); surface(&mut rng, &hit, &mut phot, &mut env, &mut data); - travel(&mut data, &mut phot, &env, &index, bump_dist); + travel(&mut data, &mut phot, &env, bump_dist); }, Event::Boundary(boundary_hit) => { - travel(&mut data, &mut phot, &env, &index, boundary_hit.dist()); + travel(&mut data, &mut phot, &env, boundary_hit.dist()); input.bound.apply(rng, &boundary_hit, &mut phot); } } diff --git a/src/sim/engine/engines/standard.rs b/src/sim/engine/engines/standard.rs index 2aaee4b..12737c8 100644 --- a/src/sim/engine/engines/standard.rs +++ b/src/sim/engine/engines/standard.rs @@ -1,9 +1,7 @@ //! Standard photon-lifetime engine function. use crate::{ - geom::Trace, - phys::Photon, - sim::{scatter::scatter, surface::surface, travel::travel, Event, Input, Output}, + geom::Trace, io::output::{self, Output}, phys::Photon, sim::{scatter::scatter, surface::surface, travel::travel, Event, Input} }; use rand::{rngs::ThreadRng, Rng}; @@ -11,11 +9,18 @@ use rand::{rngs::ThreadRng, Rng}; #[allow(clippy::expect_used)] #[inline] pub fn standard(input: &Input, mut data: &mut Output, mut rng: &mut ThreadRng, mut phot: Photon) { - // Check photon is within the grid. - if let Some(index) = input.grid.gen_index(phot.ray().pos()) { - data.emission[index] += phot.power() * phot.weight(); - } else { - panic!("Photon was not emitted within the grid."); + // Add the emission to output volumes. + // if let Some(index) = input.grid.gen_index(phot.ray().pos()) { + // data.emission[index] += phot.power() * phot.weight(); + // } else { + // panic!("Photon was not emitted within the grid."); + // } + + // Add to the emission variables in which the photon is present. + for vol in data.get_volumes_for_param_mut(output::OutputParameter::Emission) { + if let Some(index) = vol.gen_index(phot.ray().pos()) { + vol.data_mut()[index] += phot.power() * phot.weight(); + } } // Common constants. @@ -67,20 +72,20 @@ pub fn standard(input: &Input, mut data: &mut Output, mut rng: &mut ThreadRng, m // Event handling. match Event::new(voxel_dist, scat_dist, surf_hit, boundary_hit, bump_dist) { - Event::Voxel(dist) => travel(&mut data, &mut phot, &env, &index, dist + bump_dist), + Event::Voxel(dist) => travel(&mut data, &mut phot, &env, dist + bump_dist), Event::Scattering(dist) => { - travel(&mut data, &mut phot, &env, &index, dist); + travel(&mut data, &mut phot, &env,dist); scatter(&mut rng, &mut phot, &env); } Event::Surface(hit) => { - travel(&mut data, &mut phot, &env, &index, hit.dist()); + travel(&mut data, &mut phot, &env,hit.dist()); surface(&mut rng, &hit, &mut phot, &mut env, &mut data); - travel(&mut data, &mut phot, &env, &index, bump_dist); + travel(&mut data, &mut phot, &env,bump_dist); }, Event::Boundary(boundary_hit) => { - travel(&mut data, &mut phot, &env, &index, boundary_hit.dist()); + travel(&mut data, &mut phot, &env, boundary_hit.dist()); input.bound.apply(rng, &boundary_hit, &mut phot); - travel(&mut data, &mut phot, &env, &index, 100.0 * bump_dist); + travel(&mut data, &mut phot, &env, 100.0 * bump_dist); } } diff --git a/src/sim/run.rs b/src/sim/run.rs index 8d6cfae..08b063a 100644 --- a/src/sim/run.rs +++ b/src/sim/run.rs @@ -2,8 +2,9 @@ use crate::{ err::Error, - sim::{Engine, Input, Output}, + sim::{Engine, Input}, tools::ProgressBar, + io::output::Output, }; use rand::thread_rng; use rayon::prelude::*; @@ -17,8 +18,8 @@ use std::sync::{Arc, Mutex}; pub fn multi_thread<'a>( engine: &Engine, input: Input<'a>, - output: &Output<'a>, -) -> Result, Error> { + output: &Output, +) -> Result { let pb = ProgressBar::new("MCRT", input.sett.num_phot()); let pb = Arc::new(Mutex::new(pb)); @@ -36,7 +37,7 @@ pub fn multi_thread<'a>( let mut data = out.pop().expect("No data received."); while let Some(o) = out.pop() { - data += &o; + data += o; } Ok(data) @@ -49,9 +50,9 @@ pub fn multi_thread<'a>( fn thread<'a>( engine: &Engine, input: Input<'a>, - mut output: Output<'a>, + mut output: Output, pb: &Arc>, -) -> Output<'a> { +) -> Output { let mut rng = thread_rng(); let phot_energy = input.light.power() / input.sett.num_phot() as f64; diff --git a/src/sim/surface.rs b/src/sim/surface.rs index 16ab0e4..d176435 100644 --- a/src/sim/surface.rs +++ b/src/sim/surface.rs @@ -5,7 +5,8 @@ use crate::{ img::Colour, ord::cartesian::{X, Y}, phys::{Crossing, Local, Photon}, - sim::{Attribute, Output}, + sim::Attribute, + io::output::Output, }; use rand::{rngs::ThreadRng, Rng}; diff --git a/src/sim/travel.rs b/src/sim/travel.rs index ffc45f9..08100f4 100644 --- a/src/sim/travel.rs +++ b/src/sim/travel.rs @@ -2,24 +2,35 @@ use crate::{ phys::{Local, Photon}, - sim::Output, - geom::Cube, + io::output::{Output, OutputParameter}, }; use physical_constants::SPEED_OF_LIGHT_IN_VACUUM; /// Move the photon forward and record the flight. #[inline] -pub fn travel(data: &mut Output, phot: &mut Photon, env: &Local, index: &Option<([usize; 3], Cube)>, dist: f64) { +pub fn travel(data: &mut Output, phot: &mut Photon, env: &Local, dist: f64) { debug_assert!(dist > 0.0); - match *index { - Some((idx, _)) => { - let weight_power_dist = phot.weight() * phot.power() * dist; - data.energy[idx] += weight_power_dist * env.ref_index() / SPEED_OF_LIGHT_IN_VACUUM; - data.absorptions[idx] += weight_power_dist * env.abs_coeff(); - data.shifts[idx] += weight_power_dist * env.shift_coeff(); - }, - None => {}, + // Energy Density. + let weight_power_dist = phot.weight() * phot.power() * dist; + for vol in data.get_volumes_for_param_mut(OutputParameter::Energy) { + if let Some(index) = vol.gen_index(phot.ray().pos()) { + vol.data_mut()[index] += weight_power_dist * env.ref_index() / SPEED_OF_LIGHT_IN_VACUUM; + } + } + + // Absorption. + for vol in data.get_volumes_for_param_mut(OutputParameter::Absorption) { + if let Some(index) = vol.gen_index(phot.ray().pos()) { + vol.data_mut()[index] += weight_power_dist * env.abs_coeff(); + } + } + + // Shifts. + for vol in data.get_volumes_for_param_mut(OutputParameter::Shift) { + if let Some(index) = vol.gen_index(phot.ray().pos()) { + vol.data_mut()[index] += weight_power_dist * env.shift_coeff(); + } } phot.ray_mut().travel(dist); From 33754c9357f9726a0ad5871581e12be271e79221 Mon Sep 17 00:00:00 2001 From: Sam Morrell Date: Wed, 14 Aug 2024 23:34:47 +0100 Subject: [PATCH 12/16] Fixed large grid unit tests. These grids were too large and wouldn't compile on a computer with less RAM. I hadn't forgotten to change this back. --- src/io/output/output_config.rs | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/io/output/output_config.rs b/src/io/output/output_config.rs index d5a95fa..9bdee74 100644 --- a/src/io/output/output_config.rs +++ b/src/io/output/output_config.rs @@ -258,7 +258,7 @@ mod tests { let conf_str = r#" { volumes: { - full_vol: { boundary: [[0, 0, 0], [10, 10, 10]], res: [8000, 8000, 1000], param: "energy" }, + full_vol: { boundary: [[0, 0, 0], [10, 10, 10]], res: [100, 100, 100], param: "energy" }, partial_vol: { boundary: [[2.5, 2.5, 0], [2.5, 2.5, 10]], res: [100, 100, 10], param: "energy" }, }, planes: { From 242e7d82c8663c1a0876d5ccb03f5b0ebc0dfd59 Mon Sep 17 00:00:00 2001 From: Sam Morrell Date: Thu, 15 Aug 2024 02:15:51 +0100 Subject: [PATCH 13/16] First working version with implementation tests. - Added saving to new `Output` struct. - Added saving to output plane struct. - Converted voxel distance measurement in engines to new system. - Fixed bug where killed photons would continue to propagate in all engines. --- src/io/output/output.rs | 34 +++++++++++++++++++++++++- src/io/output/output_plane.rs | 8 ++++++ src/sim/engine/engines/fluorescence.rs | 11 ++++----- src/sim/engine/engines/photo.rs | 12 ++++----- src/sim/engine/engines/raman.rs | 12 ++++----- src/sim/engine/engines/standard.rs | 13 ++++------ 6 files changed, 61 insertions(+), 29 deletions(-) diff --git a/src/io/output/output.rs b/src/io/output/output.rs index 1c60226..af4b61a 100644 --- a/src/io/output/output.rs +++ b/src/io/output/output.rs @@ -103,6 +103,38 @@ impl AddAssign for Output { impl Save for Output { #[inline] fn save_data(&self, out_dir: &Path) -> Result<(), Error> { - todo!() + + for (vol, name) in self.vol.iter().zip(self.reg.vol_reg.names_list()) { + let path = out_dir.join(format!("volume_{}.nc", name.to_string())); + vol.save(&path)?; + } + + for (plane, name) in self.plane.iter().zip(self.reg.plane_reg.names_list()) { + let path = out_dir.join(format!("plane_{}.nc", name.to_string())); + plane.save(&path)?; + } + + for (name, index) in self.reg.spec_reg.set().map().iter() { + self.specs[*index].save(&out_dir.join(&format!("spectrometer_{}.csv", name)))?; + } + + for (name, index) in self.reg.img_reg.set().map().iter() { + self.imgs[*index].save(&out_dir.join(&format!("img_{}.png", name)))?; + } + + for (name, index) in self.reg.ccd_reg.set().map().iter() { + self.ccds[*index].save(&out_dir.join(&format!("ccd_{}.nc", name)))?; + } + + for (n, photo) in self.photos.iter().enumerate() { + photo.save(&out_dir.join(&format!("photo_{:03}.png", n)))?; + } + + for (name, index) in self.reg.phot_cols_reg.set().map().iter() { + self.phot_cols[*index] + .save(&out_dir.join(&format!("photon_collector_{}.csv", name)))?; + } + + Ok(()) } } diff --git a/src/io/output/output_plane.rs b/src/io/output/output_plane.rs index 8d16fc2..dcde5eb 100644 --- a/src/io/output/output_plane.rs +++ b/src/io/output/output_plane.rs @@ -3,6 +3,7 @@ use std::ops::AddAssign; use ndarray::Array2; use serde::{Deserialize, Serialize}; use crate::{ + fs::Save, math::{Point2, Point3}, ord::cartesian::{X, Y} }; @@ -106,6 +107,13 @@ impl AddAssign<&Self> for OutputPlane { } } +impl Save for OutputPlane { + fn save_data(&self, path: &std::path::Path) -> Result<(), crate::err::Error> { + self.data.save(&path)?; + Ok(()) + } +} + #[cfg(test)] mod tests { use super::*; diff --git a/src/sim/engine/engines/fluorescence.rs b/src/sim/engine/engines/fluorescence.rs index a976aa1..a25c1fb 100644 --- a/src/sim/engine/engines/fluorescence.rs +++ b/src/sim/engine/engines/fluorescence.rs @@ -71,12 +71,7 @@ pub fn fluorescence( ); // Interaction distances. - let voxel_dist = match &index { - Some((_index, voxel)) => { - voxel.dist(phot.ray()).expect("Could not determine voxel distance.") - }, - None => f64::INFINITY, - }; + let voxel_dist = data.voxel_dist(&phot); let scat_dist = -(rng.gen::()).ln() / env.inter_coeff(); let surf_hit = input .tree @@ -98,6 +93,10 @@ pub fn fluorescence( Event::Boundary(boundary_hit) => { travel(&mut data, &mut phot, &env, boundary_hit.dist()); input.bound.apply(rng, &boundary_hit, &mut phot); + // Allow for the possibility that the photon got killed at the boundary - hence don't evolve. + if phot.weight() > 0.0 { + travel(&mut data, &mut phot, &env, 100.0 * bump_dist); + } } } diff --git a/src/sim/engine/engines/photo.rs b/src/sim/engine/engines/photo.rs index 2eb9368..6a4405b 100644 --- a/src/sim/engine/engines/photo.rs +++ b/src/sim/engine/engines/photo.rs @@ -61,13 +61,7 @@ pub fn photo( } // Interaction distances. - let index = input.grid.gen_index_voxel(phot.ray().pos()); - let voxel_dist = match &index { - Some((_index, voxel)) => { - voxel.dist(phot.ray()).expect("Could not determine voxel distance.") - }, - None => f64::INFINITY, - }; + let voxel_dist = data.voxel_dist(&phot); let scat_dist = -(rng.gen::()).ln() / env.inter_coeff(); let surf_hit = input .tree @@ -105,6 +99,10 @@ pub fn photo( Event::Boundary(boundary_hit) => { travel(&mut data, &mut phot, &env, boundary_hit.dist()); input.bound.apply(rng, &boundary_hit, &mut phot); + // Allow for the possibility that the photon got killed at the boundary - hence don't evolve. + if phot.weight() > 0.0 { + travel(&mut data, &mut phot, &env, 100.0 * bump_dist); + } } } diff --git a/src/sim/engine/engines/raman.rs b/src/sim/engine/engines/raman.rs index 48244fd..70a028f 100644 --- a/src/sim/engine/engines/raman.rs +++ b/src/sim/engine/engines/raman.rs @@ -58,13 +58,7 @@ pub fn raman( } // Interaction distances. - let index = input.grid.gen_index_voxel(phot.ray().pos()); - let voxel_dist = match &index { - Some((_index, voxel)) => { - voxel.dist(phot.ray()).expect("Could not determine voxel distance.") - }, - None => f64::INFINITY, - }; + let voxel_dist = data.voxel_dist(&phot); let scat_dist = -(rng.gen::()).ln() / env.inter_coeff(); let surf_hit = input .tree @@ -94,6 +88,10 @@ pub fn raman( Event::Boundary(boundary_hit) => { travel(&mut data, &mut phot, &env, boundary_hit.dist()); input.bound.apply(rng, &boundary_hit, &mut phot); + // Allow for the possibility that the photon got killed at the boundary - hence don't evolve. + if phot.weight() > 0.0 { + travel(&mut data, &mut phot, &env, 100.0 * bump_dist); + } } } diff --git a/src/sim/engine/engines/standard.rs b/src/sim/engine/engines/standard.rs index 12737c8..a718a34 100644 --- a/src/sim/engine/engines/standard.rs +++ b/src/sim/engine/engines/standard.rs @@ -57,13 +57,7 @@ pub fn standard(input: &Input, mut data: &mut Output, mut rng: &mut ThreadRng, m } // Interaction distances. - let index = input.grid.gen_index_voxel(phot.ray().pos()); - let voxel_dist = match &index { - Some((_index, voxel)) => { - voxel.dist(phot.ray()).expect("Could not determine voxel distance.") - }, - None => f64::INFINITY, - }; + let voxel_dist = data.voxel_dist(&phot); let scat_dist = -(rng.gen::()).ln() / env.inter_coeff(); let surf_hit = input .tree @@ -85,7 +79,10 @@ pub fn standard(input: &Input, mut data: &mut Output, mut rng: &mut ThreadRng, m Event::Boundary(boundary_hit) => { travel(&mut data, &mut phot, &env, boundary_hit.dist()); input.bound.apply(rng, &boundary_hit, &mut phot); - travel(&mut data, &mut phot, &env, 100.0 * bump_dist); + // Allow for the possibility that the photon got killed at the boundary - hence don't evolve. + if phot.weight() > 0.0 { + travel(&mut data, &mut phot, &env, 100.0 * bump_dist); + } } } From 7c9bcc283c19607c6ea55973d80daa6380ce24fc Mon Sep 17 00:00:00 2001 From: Sam Morrell Date: Thu, 15 Aug 2024 17:32:27 +0100 Subject: [PATCH 14/16] Added `Display` impl for `Output` and `OutputConfig`. - Added Display implementations for supporting structs. - Fixed some warnings, mainly due to unneeded imports. --- src/data/histogram_builder.rs | 17 +++- src/img/image_builder.rs | 13 ++- src/io/output/ccd_builder.rs | 14 +++- src/io/output/output.rs | 24 ++++++ src/io/output/output_config.rs | 99 ++++++++++++++++------- src/io/output/output_plane.rs | 14 ++++ src/io/output/output_plane_builder.rs | 20 +++++ src/io/output/output_registry.rs | 2 +- src/io/output/output_volume.rs | 12 +++ src/io/output/output_volume_builder.rs | 20 +++++ src/io/output/photon_collector_builder.rs | 16 +++- src/sim/engine/engines/fluorescence.rs | 3 +- src/sim/engine/engines/photo.rs | 3 +- src/sim/engine/engines/raman.rs | 3 +- src/sim/engine/engines/standard.rs | 12 +-- 15 files changed, 224 insertions(+), 48 deletions(-) diff --git a/src/data/histogram_builder.rs b/src/data/histogram_builder.rs index 81abbb5..679b137 100644 --- a/src/data/histogram_builder.rs +++ b/src/data/histogram_builder.rs @@ -1,6 +1,9 @@ use serde::{Deserialize, Serialize}; - -use super::Histogram; +use std::fmt::{Display, Formatter}; +use crate::{ + fmt_report, + data::Histogram, +}; #[derive(Debug, Serialize, Deserialize)] @@ -14,4 +17,14 @@ impl HistogramBuilder { pub fn build(&self) -> Histogram { Histogram::new(self.min, self.max, self.bins) } +} + +impl Display for HistogramBuilder { + #[inline] + fn fmt(&self, fmt: &mut Formatter) -> Result<(), std::fmt::Error> { + fmt_report!(fmt, self.min, "min"); + fmt_report!(fmt, self.min, "max"); + fmt_report!(fmt, self.bins, "bins"); + Ok(()) + } } \ No newline at end of file diff --git a/src/img/image_builder.rs b/src/img/image_builder.rs index df33038..1baf15a 100644 --- a/src/img/image_builder.rs +++ b/src/img/image_builder.rs @@ -1,5 +1,9 @@ +use std::fmt::{Display, Formatter}; use serde::{Serialize, Deserialize}; -use crate::img::{Image, Colour}; +use crate::{ + fmt_report, + img::{Image, Colour} +}; #[derive(Debug, Serialize, Deserialize)] pub struct ImageBuilder { @@ -11,4 +15,11 @@ impl ImageBuilder { let background = Colour::new(0.0, 0.0, 0.0, 1.0); Image::new_blank(self.res, background) } +} + +impl Display for ImageBuilder { + fn fmt(&self, fmt: &mut Formatter<'_>) -> std::fmt::Result { + fmt_report!(fmt, format!("{} x {}", self.res[0], self.res[1]), "resolution"); + Ok(()) + } } \ No newline at end of file diff --git a/src/io/output/ccd_builder.rs b/src/io/output/ccd_builder.rs index a378ba7..d8d0bd4 100644 --- a/src/io/output/ccd_builder.rs +++ b/src/io/output/ccd_builder.rs @@ -1,6 +1,10 @@ +use std::fmt::{Display, Formatter}; use serde::{Serialize, Deserialize}; use ndarray::Array3; -use crate::ord::cartesian::{X, Y}; +use crate::{ + fmt_report, + ord::cartesian::{X, Y} +}; #[derive(Debug, Serialize, Deserialize)] pub struct CcdBuilder { @@ -12,4 +16,12 @@ impl CcdBuilder { pub fn build(&self) -> Array3 { Array3::zeros([self.res[X], self.res[Y], self.bins]) } +} + +impl Display for CcdBuilder { + fn fmt(&self, fmt: &mut Formatter<'_>) -> std::fmt::Result { + fmt_report!(fmt, format!("{} x {}", self.res[0], self.res[1]), "resolution"); + fmt_report!(fmt, self.bins, "bins"); + Ok(()) + } } \ No newline at end of file diff --git a/src/io/output/output.rs b/src/io/output/output.rs index af4b61a..9dcfdc5 100644 --- a/src/io/output/output.rs +++ b/src/io/output/output.rs @@ -1,4 +1,5 @@ use crate::{ + fmt_report, fs::Save, data::Histogram, img::Image, @@ -8,6 +9,7 @@ use crate::{ use std::{ ops::AddAssign, path::Path, + fmt::{Display, Formatter}, }; use ndarray::Array3; @@ -138,3 +140,25 @@ impl Save for Output { Ok(()) } } + +impl Display for Output { + #[inline] + fn fmt(&self, fmt: &mut Formatter) -> Result<(), std::fmt::Error> { + writeln!(fmt, "...")?; + + fmt_report!(fmt, self.reg.vol_reg, "output volume register"); + fmt_report!(fmt, self.reg.plane_reg, "output plane register"); + fmt_report!(fmt, self.reg.phot_cols_reg, "photon collector register"); + fmt_report!(fmt, self.reg.spec_reg, "spectrometer register"); + fmt_report!(fmt, self.reg.img_reg, "imager register"); + fmt_report!(fmt, self.reg.ccd_reg, "ccd register"); + + fmt_report!(fmt, self.specs.len(), "spectrometers"); + fmt_report!(fmt, self.imgs.len(), "images"); + fmt_report!(fmt, self.ccds.len(), "ccds"); + + fmt_report!(fmt, self.photos.len(), "photos"); + fmt_report!(fmt, self.phot_cols.len(), "photon collectors"); + Ok(()) + } +} \ No newline at end of file diff --git a/src/io/output/output_config.rs b/src/io/output/output_config.rs index 9bdee74..5e45395 100644 --- a/src/io/output/output_config.rs +++ b/src/io/output/output_config.rs @@ -6,6 +6,7 @@ use std::{ use serde::Serialize; use arctk_attr::file; use crate::{ + fmt_report, data::HistogramBuilder, img::ImageBuilder, io::output::{OutputPlaneBuilder, OutputVolumeBuilder, PhotonCollectorBuilder, Output}, @@ -164,39 +165,81 @@ impl OutputConfig { } impl fmt::Display for OutputConfig { - fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result { - write!(f, "OutputConfig {{\n")?; - - if let Some(volumes) = &self.volumes { - write!(f, " volumes: {{\n")?; - for (name, volume) in volumes { - write!(f, " {}: {:?}\n", name, volume)?; - } - write!(f, " }}\n")?; + #[inline] + fn fmt(&self, fmt: &mut fmt::Formatter) -> Result<(), std::fmt::Error> { + writeln!(fmt, "...")?; + + match &self.volumes { + Some(vols) => { + fmt_report!(fmt, "...", "volume outputs"); + for (key, vol) in vols { + fmt_report!(fmt, vol, key); + } + }, + None => fmt_report!(fmt, "none", "volume outputs") } - - if let Some(planes) = &self.planes { - write!(f, " planes: {{\n")?; - for (name, plane) in planes { - write!(f, " {}: {:?}\n", name, plane)?; - } - write!(f, " }}\n")?; + + match &self.planes { + Some(planes) => { + fmt_report!(fmt, "...", "plane outputs"); + for (key, plane) in planes { + fmt_report!(fmt, plane, key); + } + }, + None => fmt_report!(fmt, "none", "plane outputs") } - - if let Some(photon_collectors) = &self.photon_collectors { - write!(f, " photon_collectors: {{\n")?; - for (name, collector) in photon_collectors { - write!(f, " {}: {:?}\n", name, collector)?; - } - write!(f, " }}\n")?; + + match &self.photon_collectors { + Some(pcs) => { + fmt_report!(fmt, "...", "photon collectors"); + for (key, pc) in pcs { + fmt_report!(fmt, pc, key); + } + }, + None => fmt_report!(fmt, "none", "photon collectors") } - // TODO: Add output for spectra. - // TODO: Add output for images. - // TODO: Add output for CCDs. - // TODO: Add output for photos. + match &self.spectra { + Some(spectra) => { + fmt_report!(fmt, "...", "spectra"); + for (key, spec) in spectra { + fmt_report!(fmt, spec, key); + } + }, + None => fmt_report!(fmt, "none", "spectra") + } + + match &self.images { + Some(imgs) => { + fmt_report!(fmt, "...", "images"); + for (key, img) in imgs { + fmt_report!(fmt, img, key); + } + }, + None => fmt_report!(fmt, "none", "images") + } + + match &self.ccds { + Some(ccds) => { + fmt_report!(fmt, "...", "ccds"); + for (key, ccd) in ccds { + fmt_report!(fmt, ccd, key); + } + }, + None => fmt_report!(fmt, "none", "ccds") + } + + match &self.photos { + Some(photos) => { + fmt_report!(fmt, "...", "photos"); + for (key, photo) in photos { + fmt_report!(fmt, photo, key); + } + }, + None => fmt_report!(fmt, "none", "photos") + } - write!(f, "}}") + Ok(()) } } diff --git a/src/io/output/output_plane.rs b/src/io/output/output_plane.rs index dcde5eb..c2153a3 100644 --- a/src/io/output/output_plane.rs +++ b/src/io/output/output_plane.rs @@ -2,6 +2,7 @@ use std::ops::AddAssign; use ndarray::Array2; use serde::{Deserialize, Serialize}; +use std::fmt; use crate::{ fs::Save, math::{Point2, Point3}, @@ -17,6 +18,17 @@ pub enum AxisAlignedPlane { YZ, } +impl fmt::Display for AxisAlignedPlane { + fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result { + let plane_str = match self { + AxisAlignedPlane::XY => "XY", + AxisAlignedPlane::XZ => "XZ", + AxisAlignedPlane::YZ => "YZ", + }; + write!(f, "{}", plane_str) + } +} + impl AxisAlignedPlane { pub fn project_onto_plane(&self, p: &Point3) -> (f64, f64) { match *self { @@ -27,6 +39,8 @@ impl AxisAlignedPlane { } } + + #[derive(Debug, Clone)] pub struct OutputPlane { mins: Point2, diff --git a/src/io/output/output_plane_builder.rs b/src/io/output/output_plane_builder.rs index d2ef79a..fbe5eba 100644 --- a/src/io/output/output_plane_builder.rs +++ b/src/io/output/output_plane_builder.rs @@ -1,7 +1,10 @@ +use std::fmt::{Display, Formatter}; use serde::{Deserialize, Serialize}; use crate::{ + fmt_report, math::Point2, io::output::{OutputPlane, AxisAlignedPlane}, + ord::cartesian::{X, Y}, }; #[derive(Debug, Deserialize, Serialize)] @@ -36,3 +39,20 @@ mod tests { assert_eq!(outvol.pix_area(), 1.0); } } + +impl Display for OutputPlaneBuilder { + #[inline] + fn fmt(&self, fmt: &mut Formatter) -> Result<(), std::fmt::Error> { + writeln!(fmt, "...")?; + fmt_report!(fmt, "...", "boundary"); + fmt_report!(fmt, format!("[{}, {}", self.boundary.0.x(), self.boundary.0.y()), "mins"); + fmt_report!(fmt, format!("[{}, {}", self.boundary.0.x(), self.boundary.0.y()), "maxs"); + fmt_report!( + fmt, + &format!("[{} x {}]", self.res[X], self.res[Y]), + "resolution" + ); + fmt_report!(fmt, self.plane, "parameter"); + Ok(()) + } +} \ No newline at end of file diff --git a/src/io/output/output_registry.rs b/src/io/output/output_registry.rs index 6e0a57b..760d1a3 100644 --- a/src/io/output/output_registry.rs +++ b/src/io/output/output_registry.rs @@ -1,5 +1,5 @@ use crate::{ - io::output::OutputConfig, ord::{Name, Register} + io::output::OutputConfig, ord::Register, }; #[derive(Clone)] diff --git a/src/io/output/output_volume.rs b/src/io/output/output_volume.rs index ca5fed5..8394152 100644 --- a/src/io/output/output_volume.rs +++ b/src/io/output/output_volume.rs @@ -25,6 +25,18 @@ pub enum OutputParameter { Shift, } +impl std::fmt::Display for OutputParameter { + fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result { + let param_str = match self { + OutputParameter::Emission => "Emission", + OutputParameter::Energy => "Energy", + OutputParameter::Absorption => "Absorption", + OutputParameter::Shift => "Shift", + }; + write!(f, "{}", param_str) + } +} + #[derive(Debug, Clone)] pub struct OutputVolume { boundary: Cube, diff --git a/src/io/output/output_volume_builder.rs b/src/io/output/output_volume_builder.rs index af6b234..f926ae4 100644 --- a/src/io/output/output_volume_builder.rs +++ b/src/io/output/output_volume_builder.rs @@ -1,8 +1,11 @@ use serde::{Deserialize, Serialize}; +use std::fmt::{Display, Formatter}; use crate::{ + fmt_report, geom::Cube, math::Vec3, io::output::{OutputVolume, OutputParameter}, + ord::cartesian::{X, Y, Z}, }; /// Configuration for the OutputVolume. @@ -22,6 +25,23 @@ impl OutputVolumeBuilder { } } +impl Display for OutputVolumeBuilder { + #[inline] + fn fmt(&self, fmt: &mut Formatter) -> Result<(), std::fmt::Error> { + writeln!(fmt, "...")?; + fmt_report!(fmt, "...", "boundary"); + fmt_report!(fmt, format!("[{}, {}", self.boundary.0.x(), self.boundary.0.y()), "mins"); + fmt_report!(fmt, format!("[{}, {}", self.boundary.0.x(), self.boundary.0.y()), "maxs"); + fmt_report!( + fmt, + &format!("[{} x {} x {}]", self.res[X], self.res[Y], self.res[Z]), + "resolution" + ); + fmt_report!(fmt, self.param, "parameter"); + Ok(()) + } +} + #[cfg(test)] mod tests { use super::*; diff --git a/src/io/output/photon_collector_builder.rs b/src/io/output/photon_collector_builder.rs index b5dc20a..c4d18be 100644 --- a/src/io/output/photon_collector_builder.rs +++ b/src/io/output/photon_collector_builder.rs @@ -1,5 +1,6 @@ +use std::fmt::{Display, Formatter}; use serde::{Serialize, Deserialize}; -use crate::io::output::PhotonCollector; +use crate::{fmt_report, io::output::PhotonCollector}; #[derive(Default, Debug, Serialize, Deserialize)] pub struct PhotonCollectorBuilder { @@ -17,6 +18,19 @@ impl PhotonCollectorBuilder { } } +impl Display for PhotonCollectorBuilder { + #[inline] + fn fmt(&self, fmt: &mut Formatter<'_>) -> std::fmt::Result { + let kill_str = match self.kill_photons { + Some(kf) => kf, + None => false, + }; + writeln!(fmt, "PhotonCollector: ")?; + fmt_report!(fmt, kill_str, "kill on collect"); + Ok(()) + } +} + #[cfg(test)] mod tests { use super::*; diff --git a/src/sim/engine/engines/fluorescence.rs b/src/sim/engine/engines/fluorescence.rs index a25c1fb..b9bd9d0 100644 --- a/src/sim/engine/engines/fluorescence.rs +++ b/src/sim/engine/engines/fluorescence.rs @@ -1,7 +1,6 @@ //! Fluorescence photon-lifetime engine function. use crate::{ - geom::Trace, math::Formula, phys::{Local, Photon}, sim::{scatter::scatter, surface::surface, travel::travel, Event, Input}, @@ -95,7 +94,7 @@ pub fn fluorescence( input.bound.apply(rng, &boundary_hit, &mut phot); // Allow for the possibility that the photon got killed at the boundary - hence don't evolve. if phot.weight() > 0.0 { - travel(&mut data, &mut phot, &env, 100.0 * bump_dist); + travel(&mut data, &mut phot, &env, bump_dist); } } } diff --git a/src/sim/engine/engines/photo.rs b/src/sim/engine/engines/photo.rs index 6a4405b..24d3e74 100644 --- a/src/sim/engine/engines/photo.rs +++ b/src/sim/engine/engines/photo.rs @@ -1,7 +1,6 @@ //! Photography photon-lifetime engine function. use crate::{ - geom::Trace, img::Colour, phys::Photon, io::output::{Output, OutputParameter}, @@ -101,7 +100,7 @@ pub fn photo( input.bound.apply(rng, &boundary_hit, &mut phot); // Allow for the possibility that the photon got killed at the boundary - hence don't evolve. if phot.weight() > 0.0 { - travel(&mut data, &mut phot, &env, 100.0 * bump_dist); + travel(&mut data, &mut phot, &env, bump_dist); } } } diff --git a/src/sim/engine/engines/raman.rs b/src/sim/engine/engines/raman.rs index 70a028f..02d2ee8 100644 --- a/src/sim/engine/engines/raman.rs +++ b/src/sim/engine/engines/raman.rs @@ -1,7 +1,6 @@ //! Raman photon-lifetime engine function. use crate::{ - geom::Trace, math::Point3, phys::Photon, sim::{scatter::shift_scatter, surface::surface, travel::travel, Event, Input}, @@ -90,7 +89,7 @@ pub fn raman( input.bound.apply(rng, &boundary_hit, &mut phot); // Allow for the possibility that the photon got killed at the boundary - hence don't evolve. if phot.weight() > 0.0 { - travel(&mut data, &mut phot, &env, 100.0 * bump_dist); + travel(&mut data, &mut phot, &env, bump_dist); } } } diff --git a/src/sim/engine/engines/standard.rs b/src/sim/engine/engines/standard.rs index a718a34..09fd542 100644 --- a/src/sim/engine/engines/standard.rs +++ b/src/sim/engine/engines/standard.rs @@ -1,7 +1,9 @@ //! Standard photon-lifetime engine function. use crate::{ - geom::Trace, io::output::{self, Output}, phys::Photon, sim::{scatter::scatter, surface::surface, travel::travel, Event, Input} + io::output::{self, Output}, + phys::Photon, + sim::{scatter::scatter, surface::surface, travel::travel, Event, Input}, }; use rand::{rngs::ThreadRng, Rng}; @@ -9,12 +11,6 @@ use rand::{rngs::ThreadRng, Rng}; #[allow(clippy::expect_used)] #[inline] pub fn standard(input: &Input, mut data: &mut Output, mut rng: &mut ThreadRng, mut phot: Photon) { - // Add the emission to output volumes. - // if let Some(index) = input.grid.gen_index(phot.ray().pos()) { - // data.emission[index] += phot.power() * phot.weight(); - // } else { - // panic!("Photon was not emitted within the grid."); - // } // Add to the emission variables in which the photon is present. for vol in data.get_volumes_for_param_mut(output::OutputParameter::Emission) { @@ -81,7 +77,7 @@ pub fn standard(input: &Input, mut data: &mut Output, mut rng: &mut ThreadRng, m input.bound.apply(rng, &boundary_hit, &mut phot); // Allow for the possibility that the photon got killed at the boundary - hence don't evolve. if phot.weight() > 0.0 { - travel(&mut data, &mut phot, &env, 100.0 * bump_dist); + travel(&mut data, &mut phot, &env, bump_dist); } } } From b8714410867e74665edd5229ab01254ce69c3678 Mon Sep 17 00:00:00 2001 From: Sam Morrell Date: Thu, 15 Aug 2024 22:55:35 +0100 Subject: [PATCH 15/16] Added boundary conditions config to MCRT. Removed the grid config and replaced with the boundary conditions. Many changes to support these changes. Added structs and types to support deserialisation and building. Integrated into the MCRT binary. Note, this commit breaks the fluorescence engine, but I think we can work around this. --- src/bin/mcrt.rs | 72 ++++----- src/geom/domain/boundary.rs | 17 +- src/geom/domain/boundary_builder.rs | 153 ++++++++++++++++++ src/geom/domain/mod.rs | 3 +- src/phys/mod.rs | 3 +- src/phys/reflectance.rs | 4 +- src/phys/reflectance_builder.rs | 54 +++++++ ...bute_linker_linker_linker_linker_linker.rs | 48 ++---- src/sim/engine/engines/fluorescence.rs | 21 ++- src/sim/input.rs | 7 +- src/sim/param/parameters.rs | 12 +- src/sim/param/parameters_builder.rs | 16 +- src/sim/param/parameters_builder_loader.rs | 10 +- 13 files changed, 301 insertions(+), 119 deletions(-) create mode 100644 src/geom/domain/boundary_builder.rs create mode 100644 src/phys/reflectance_builder.rs diff --git a/src/bin/mcrt.rs b/src/bin/mcrt.rs index 88a38df..25729f9 100644 --- a/src/bin/mcrt.rs +++ b/src/bin/mcrt.rs @@ -8,12 +8,11 @@ use std::{ use aetherus::{ args, fs::{File, Load, Save}, - geom::{Tree, Boundary}, - ord::{Build, Link, Register, Set}, + geom::Tree, + ord::{Build, Link}, report, - io::output::{Output, PhotonCollector}, sim::{ - run, AttributeLinkerLinkerLinkerLinkerLinker as Attr, Engine, Input, Parameters, + run, Input, Parameters, ParametersBuilderLoader, }, util::{ @@ -22,7 +21,6 @@ use aetherus::{ fmt::term, }, }; -use nalgebra::base; /// Backup print width if the terminal width can not be determined. const BACKUP_TERM_WIDTH: usize = 80; @@ -41,13 +39,13 @@ fn main() { report!(engine, "engine"); let sett = params.sett; report!(sett, "settings"); - let grid = params.grid; - report!(grid, "measurement grid"); + let bound = params.boundary; + report!(bound, "boundary"); let mats = params.mats; report!(mats, "materials"); - sub_section(term_width, "Registration"); - let (spec_reg, img_reg, ccd_reg, phot_col_reg) = gen_detector_registers(¶ms.attrs); + //sub_section(term_width, "Registration"); + //let (spec_reg, img_reg, ccd_reg, phot_col_reg) = gen_detector_registers(¶ms.attrs); // let base_output = gen_base_output( // &engine, // &grid, @@ -90,8 +88,6 @@ fn main() { * For now we hard-code this to kill, but we can link this to configuration soon. * TODO: We probably want to implement the MPI adjacent rank transfer here too. */ - let bound = Boundary::new_kill(grid.boundary().clone()); - report!(bound, "boundary conditions"); sub_section(term_width, "Growing"); let tree = Tree::new(¶ms.tree, &surfs); @@ -104,7 +100,7 @@ fn main() { .fold(base_output.clone(), |mut output, (light_idx, (light_id, light))| { section(term_width, &format!("Running for light {} ({} / {})", light_id, light_idx + 1, nlights)); report!(light, light_id); - let input = Input::new(&spec_reg, &mats, &attrs, light, &tree, &grid, &sett, &bound); + let input = Input::new(&base_output.reg.spec_reg, &mats, &attrs, light, &tree, &sett, &bound); let data = run::multi_thread(&engine, input, &base_output).expect("Failed to run MCRT."); @@ -180,37 +176,37 @@ fn load_parameters(term_width: usize, in_dir: &Path, params_path: &Path) -> Para params } -/// Generate the detector registers. -fn gen_detector_registers(attrs: &Set) -> (Register, Register, Register, Register) { - let mut spec_names = Vec::new(); - let mut img_names = Vec::new(); - let mut ccd_names = Vec::new(); - let mut phot_col_names = Vec::new(); - - for attr in attrs.map().values() { - match *attr { - Attr::Spectrometer(ref name, ..) => spec_names.push(name.clone()), - Attr::Imager(ref name, ..) => img_names.push(name.clone()), - Attr::Ccd(ref name, ..) => ccd_names.push(name.clone()), - Attr::PhotonCollector(ref name, ..) => phot_col_names.push(name.clone()), - _ => {} - } - } +// /// Generate the detector registers. +// fn gen_detector_registers(attrs: &Set) -> (Register, Register, Register, Register) { +// let mut spec_names = Vec::new(); +// let mut img_names = Vec::new(); +// let mut ccd_names = Vec::new(); +// let mut phot_col_names = Vec::new(); + +// for attr in attrs.map().values() { +// match *attr { +// Attr::Spectrometer(ref name, ..) => spec_names.push(name.clone()), +// Attr::Imager(ref name, ..) => img_names.push(name.clone()), +// Attr::Ccd(ref name, ..) => ccd_names.push(name.clone()), +// Attr::PhotonCollector(ref name, ..) => phot_col_names.push(name.clone()), +// _ => {} +// } +// } - let spec_reg = Register::new(spec_names); - report!(spec_reg, "spectrometer register"); +// let spec_reg = Register::new(spec_names); +// report!(spec_reg, "spectrometer register"); - let img_reg = Register::new(img_names); - report!(img_reg, "imager register"); +// let img_reg = Register::new(img_names); +// report!(img_reg, "imager register"); - let ccd_reg = Register::new(ccd_names); - report!(ccd_reg, "ccd register"); +// let ccd_reg = Register::new(ccd_names); +// report!(ccd_reg, "ccd register"); - let phot_col_reg = Register::new(phot_col_names); - report!(phot_col_reg, "photon collector register"); +// let phot_col_reg = Register::new(phot_col_names); +// report!(phot_col_reg, "photon collector register"); - (spec_reg, img_reg, ccd_reg, phot_col_reg) -} +// (spec_reg, img_reg, ccd_reg, phot_col_reg) +// } // Generate the base output instance. // fn gen_base_output<'a>( diff --git a/src/geom/domain/boundary.rs b/src/geom/domain/boundary.rs index ddad49a..024ec5f 100644 --- a/src/geom/domain/boundary.rs +++ b/src/geom/domain/boundary.rs @@ -11,14 +11,15 @@ use std::fmt::{Display, Formatter}; /// Struct that represents a boundary. /// This will be used to determine how the boundary conditions behaves when it interacts /// with photon packets. +#[derive(Clone)] pub struct Boundary { - bounding_box: Cube, - top: BoundaryCondition, - bottom: BoundaryCondition, - north: BoundaryCondition, - east: BoundaryCondition, - south: BoundaryCondition, - west: BoundaryCondition, + pub bounding_box: Cube, + pub top: BoundaryCondition, + pub bottom: BoundaryCondition, + pub north: BoundaryCondition, + pub east: BoundaryCondition, + pub south: BoundaryCondition, + pub west: BoundaryCondition, } impl Boundary { @@ -299,7 +300,7 @@ impl Display for BoundaryDirection { } } -#[derive(Default, Debug, PartialEq)] +#[derive(Default, Debug, PartialEq, Clone)] pub enum BoundaryCondition { /// Any photon packet that intersects with this boundary will be down-weighted /// and removed from the simulation. diff --git a/src/geom/domain/boundary_builder.rs b/src/geom/domain/boundary_builder.rs new file mode 100644 index 0000000..d9deafb --- /dev/null +++ b/src/geom/domain/boundary_builder.rs @@ -0,0 +1,153 @@ +use serde::{Deserialize, Serialize}; +use std::fmt::{Display, Formatter}; +use arctk_attr::file; +use crate::{ + fmt_report, + geom::{Boundary, BoundaryCondition, Cube}, + math::Vec3, + phys::{ReflectanceBuilder, ReflectanceBuilderShim}, +}; + +#[file] +#[derive(Serialize)] +pub struct BoundaryBuilder { + boundary: (Vec3, Vec3), + top: Option, + bottom: Option, + north: Option, + east: Option, + south: Option, + west: Option, +} + +impl BoundaryBuilder { + pub fn build(&self) -> Boundary { + let bounding_box = Cube::new(self.boundary.0.data().into(), self.boundary.1.data().into()); + let top = match &self.top { + Some(a) => a.build(), + None => BoundaryCondition::default(), + }; + let bottom: BoundaryCondition = match &self.bottom { + Some(a) => a.build(), + None => BoundaryCondition::default(), + }; + let north: BoundaryCondition = match &self.north { + Some(a) => a.build(), + None => BoundaryCondition::default(), + }; + let south = match &self.south { + Some(a) => a.build(), + None => BoundaryCondition::default(), + }; + let east = match &self.east { + Some(a) => a.build(), + None => BoundaryCondition::default(), + }; + let west = match &self.west { + Some(a) => a.build(), + None => BoundaryCondition::default(), + }; + + Boundary { + bounding_box, + top, + bottom, + north, + south, + east, + west, + } + } +} + +impl Display for BoundaryBuilder { + fn fmt(&self, fmt: &mut Formatter<'_>) -> std::fmt::Result { + writeln!(fmt, "...")?; + fmt_report!(fmt, "...", "boundary"); + fmt_report!(fmt, format!("[{}, {}", self.boundary.0.x(), self.boundary.0.y()), "mins"); + fmt_report!(fmt, format!("[{}, {}", self.boundary.0.x(), self.boundary.0.y()), "maxs"); + + match &self.top { + Some(a) => fmt_report!(fmt, a, "top"), + None => fmt_report!(fmt, "none", "top"), + }; + + match &self.bottom { + Some(a) => fmt_report!(fmt, a, "bottom"), + None => fmt_report!(fmt, "none", "bottom"), + }; + + match &self.north { + Some(a) => fmt_report!(fmt, a, "north"), + None => fmt_report!(fmt, "none", "north"), + }; + + match &self.south { + Some(a) => fmt_report!(fmt, a, "south"), + None => fmt_report!(fmt, "none", "south"), + }; + + match &self.east { + Some(a) => fmt_report!(fmt, a, "east"), + None => fmt_report!(fmt, "none", "east"), + }; + + match &self.west { + Some(a) => fmt_report!(fmt, a, "west"), + None => fmt_report!(fmt, "none", "west"), + }; + + Ok(()) + } +} + +#[derive(Debug, Serialize, Deserialize)] +pub enum BoundaryConditionBuilder { + Kill, + Reflect(ReflectanceBuilderShim), + Periodic(f64), + #[cfg(feature = "mpi")] + MpiRank(usize), +} + +impl BoundaryConditionBuilder { + pub fn build(&self) -> BoundaryCondition { + match self { + Self::Kill => BoundaryCondition::Kill, + Self::Periodic(dist) => BoundaryCondition::Periodic(dist.clone()), + Self::Reflect(ref_shim) => { + let ref_build: ReflectanceBuilder = ref_shim.clone().into(); + let ref_model = ref_build.build().expect("Unable to load reflectance model for boundary. "); + BoundaryCondition::Reflect(ref_model) + } + } + } +} + + +impl Display for BoundaryConditionBuilder { + fn fmt(&self, fmt: &mut Formatter<'_>) -> std::fmt::Result { + match *self { + Self::Kill => { + writeln!(fmt, "Kill: ...")?; + Ok(()) + } + Self::Reflect(ref reflectance) => { + writeln!(fmt, "Reflector: ...")?; + fmt_report!(fmt, format!("{:?}, {:?}, {:?}", reflectance.0, reflectance.1, reflectance.2), "reflectance"); + Ok(()) + }, + Self::Periodic(padding) => { + writeln!(fmt, "Periodic: ...")?; + fmt_report!(fmt, padding, "padding"); + Ok(()) + }, + #[cfg(feature = "mpi")] + Self::MpiRank(rank) => { + writeln!(fmt, "MPI Rank Transfer: ...")?; + fmt_report!(fmt, padding, "destination rank"); + Ok(()) + } + } + } +} diff --git a/src/geom/domain/mod.rs b/src/geom/domain/mod.rs index 3d74bb7..4ff9404 100644 --- a/src/geom/domain/mod.rs +++ b/src/geom/domain/mod.rs @@ -1,6 +1,7 @@ //! Domain module. pub mod boundary; +pub mod boundary_builder; pub mod grid; pub mod grid_builder; pub mod surface; @@ -10,6 +11,6 @@ pub mod tree; pub mod tree_settings; pub use self::{ - boundary::*, grid::*, grid_builder::*, surface::*, surface_linker::*, + boundary::*, boundary_builder::*, grid::*, grid_builder::*, surface::*, surface_linker::*, surface_linker_loader::*, tree::*, tree_settings::*, }; diff --git a/src/phys/mod.rs b/src/phys/mod.rs index a01351d..0f12ef7 100644 --- a/src/phys/mod.rs +++ b/src/phys/mod.rs @@ -11,6 +11,7 @@ pub mod spectrum; // Builders pub mod light_linker_builder; pub mod material_builder; +pub mod reflectance_builder; pub mod spectrum_builder; // Linkers @@ -22,5 +23,5 @@ pub mod light_linker_builder_loader; pub use self::{ crossing::*, light::*, light_linker::*, light_linker_builder::*, light_linker_builder_loader::*, local::*, material::*, material_builder::*, photon::*, - reflectance::*, spectrum::*, spectrum_builder::*, + reflectance::*, reflectance_builder::*, spectrum::*, spectrum_builder::*, }; diff --git a/src/phys/reflectance.rs b/src/phys/reflectance.rs index 4e51d70..1b4f71e 100644 --- a/src/phys/reflectance.rs +++ b/src/phys/reflectance.rs @@ -2,14 +2,12 @@ use crate::{ core::Real, fmt_report, geom::{Hit, Ray}, - phys::Spectrum, + phys::{Spectrum, Photon}, sim::Attribute, }; use rand::Rng; use std::{f64::consts::PI, fmt::Display}; -use super::Photon; - /// A small utility function that checks that the provided spectrum is valid as a /// reflectance spectrum. This means that it should have values that are between 0.0 /// and 1.0. diff --git a/src/phys/reflectance_builder.rs b/src/phys/reflectance_builder.rs new file mode 100644 index 0000000..e2be5f9 --- /dev/null +++ b/src/phys/reflectance_builder.rs @@ -0,0 +1,54 @@ +use crate::{ + err::Error, + phys::SpectrumBuilder +}; + +use super::Reflectance; + +pub type ReflectanceBuilderShim = ( + Option, + Option, + Option, +); + +pub struct ReflectanceBuilder { + pub diff_ref: Option, + pub spec_ref: Option, + pub specularity: Option, +} + +impl ReflectanceBuilder { + pub fn build(&self) -> Result { + let ref_model = if self.diff_ref.is_some() { + if self.spec_ref.is_some() { + // Check that the specularity of the reflector is defined. + assert!(self.specularity.is_some()); + Reflectance::Composite { + diffuse_refspec: self.diff_ref.clone().unwrap().build()?, + specular_refspec: self.spec_ref.clone().unwrap().build()?, + specularity: self.specularity.unwrap(), + } + } else { + Reflectance::Lambertian { + refspec: self.diff_ref.clone().unwrap().build()?, + } + } + } else { + Reflectance::Specular { + refspec: self.spec_ref.clone().unwrap().build()?, + } + }; + + Ok(ref_model) + } +} + +impl From for ReflectanceBuilder { + fn from(value: ReflectanceBuilderShim) -> Self { + ReflectanceBuilder { + diff_ref: value.0, + spec_ref: value.1, + specularity: value.2, + } + } +} \ No newline at end of file diff --git a/src/sim/attribute/attribute_linker_linker_linker_linker_linker.rs b/src/sim/attribute/attribute_linker_linker_linker_linker_linker.rs index 92e1491..e303806 100644 --- a/src/sim/attribute/attribute_linker_linker_linker_linker_linker.rs +++ b/src/sim/attribute/attribute_linker_linker_linker_linker_linker.rs @@ -4,8 +4,8 @@ use crate::{ err::Error, fmt_report, math::{Point3, Vec3}, - ord::{Link, Name, Set, cartesian::{X, Y}}, - phys::{Reflectance, SpectrumBuilder}, + ord::{cartesian::{X, Y}, Link, Name, Set}, + phys::{ReflectanceBuilder, ReflectanceBuilderShim}, sim::attribute::AttributeLinkerLinkerLinkerLinker, tools::{Binner, Range}, }; @@ -28,11 +28,7 @@ pub enum AttributeLinkerLinkerLinkerLinkerLinker { Ccd(Name, [usize; 2], f64, Point3, Vec3, Binner), /// A purely reflecting material, with a provided reflectance model. /// The first coefficient is diffuse albedo, the second is specular. - Reflector( - Option, - Option, - Option, - ), + Reflector(ReflectanceBuilderShim), /// A photon collector, which collects the photon that interact with the linked entities. /// These photons can be optionally killed, or left to keep propogating. PhotonCollector(Name, bool), @@ -60,27 +56,9 @@ impl<'a> Link<'a, usize> for AttributeLinkerLinkerLinkerLinkerLinker { Self::Ccd(id, _resolution, width, center, forward, binner) => { Self::Inst::Ccd(id, _resolution, width, center, forward, binner) } - Self::Reflector(diff_ref, spec_ref, specularity) => { - let ref_model = if diff_ref.is_some() { - if spec_ref.is_some() { - // Check that the specularity of the reflector is defined. - assert!(specularity.is_some()); - Reflectance::Composite { - diffuse_refspec: diff_ref.unwrap().build()?, - specular_refspec: spec_ref.unwrap().build()?, - specularity: specularity.unwrap(), - } - } else { - Reflectance::Lambertian { - refspec: diff_ref.unwrap().build()?, - } - } - } else { - Reflectance::Specular { - refspec: spec_ref.unwrap().build()?, - } - }; - + Self::Reflector(ref_sim) => { + let ref_build: ReflectanceBuilder = ref_sim.into(); + let ref_model = ref_build.build()?; Self::Inst::Reflector(ref_model) } Self::PhotonCollector(ref id, _kill_photons) => { @@ -130,12 +108,12 @@ impl Display for AttributeLinkerLinkerLinkerLinkerLinker { fmt_report!(fmt, binner, "binner"); Ok(()) } - Self::Reflector(ref diff_ref, ref spec_ref, ref specularity) => { + Self::Reflector(ref ref_shim) => { writeln!(fmt, "Reflector: ...")?; fmt_report!( fmt, - if diff_ref.is_some() { - format!("{}", diff_ref.as_ref().unwrap()) + if ref_shim.0.is_some() { + format!("{}", ref_shim.0.as_ref().unwrap()) } else { String::from("none") }, @@ -143,8 +121,8 @@ impl Display for AttributeLinkerLinkerLinkerLinkerLinker { ); fmt_report!( fmt, - if spec_ref.is_some() { - format!("{}", spec_ref.as_ref().unwrap()) + if ref_shim.1.is_some() { + format!("{}", ref_shim.1.as_ref().unwrap()) } else { String::from("none") }, @@ -152,8 +130,8 @@ impl Display for AttributeLinkerLinkerLinkerLinkerLinker { ); fmt_report!( fmt, - if specularity.is_some() { - format!("{}", specularity.as_ref().unwrap()) + if ref_shim.2.is_some() { + format!("{}", ref_shim.2.as_ref().unwrap()) } else { String::from("none") }, diff --git a/src/sim/engine/engines/fluorescence.rs b/src/sim/engine/engines/fluorescence.rs index b9bd9d0..5864f7d 100644 --- a/src/sim/engine/engines/fluorescence.rs +++ b/src/sim/engine/engines/fluorescence.rs @@ -1,10 +1,10 @@ -//! Fluorescence photon-lifetime engine function. +// Fluorescence photon-lifetime engine function. use crate::{ + io::output::{Output, OutputParameter}, math::Formula, phys::{Local, Photon}, sim::{scatter::scatter, surface::surface, travel::travel, Event, Input}, - io::output::{Output, OutputParameter}, }; use ndarray::Array3; use rand::{rngs::ThreadRng, Rng}; @@ -20,7 +20,7 @@ pub fn fluorescence( mut rng: &mut ThreadRng, mut phot: Photon, ) { - // Add to the emission variables in which the photon is present. + // Add to the emission variables in which the photon is present. for vol in data.get_volumes_for_param_mut(OutputParameter::Emission) { if let Some(index) = vol.gen_index(phot.ray().pos()) { vol.data_mut()[index] += phot.power() * phot.weight(); @@ -60,12 +60,14 @@ pub fn fluorescence( } // Local variable modifications. - let index = input.grid.gen_index_voxel(phot.ray().pos()); + // TODO: I have had to remove this for now, as I've removed the fixed grid. + //let index = input.grid.gen_index_voxel(phot.ray().pos()); + let index = [0, 0, 0]; env = Local::new( local.ref_index(), local.scat_coeff(), local.abs_coeff(), - mu_shift.mul_add(flu_concs[index.as_ref().unwrap().0], local.shift_coeff()), + mu_shift.mul_add(flu_concs[index], local.shift_coeff()), local.asym(), ); @@ -75,7 +77,10 @@ pub fn fluorescence( let surf_hit = input .tree .scan(phot.ray().clone(), bump_dist, voxel_dist.min(scat_dist)); - let boundary_hit = input.bound.dist_boundary(phot.ray()).expect("Photon not contained in boundary. "); + let boundary_hit = input + .bound + .dist_boundary(phot.ray()) + .expect("Photon not contained in boundary. "); // Event handling. match Event::new(voxel_dist, scat_dist, surf_hit, boundary_hit, bump_dist) { @@ -88,11 +93,11 @@ pub fn fluorescence( travel(&mut data, &mut phot, &env, hit.dist()); surface(&mut rng, &hit, &mut phot, &mut local, &mut data); travel(&mut data, &mut phot, &env, bump_dist); - }, + } Event::Boundary(boundary_hit) => { travel(&mut data, &mut phot, &env, boundary_hit.dist()); input.bound.apply(rng, &boundary_hit, &mut phot); - // Allow for the possibility that the photon got killed at the boundary - hence don't evolve. + // Allow for the possibility that the photon got killed at the boundary - hence don't evolve. if phot.weight() > 0.0 { travel(&mut data, &mut phot, &env, bump_dist); } diff --git a/src/sim/input.rs b/src/sim/input.rs index 0774d0f..80af60b 100644 --- a/src/sim/input.rs +++ b/src/sim/input.rs @@ -2,7 +2,7 @@ use crate::{ fmt_report, - geom::{Boundary, Grid, Tree}, + geom::{Boundary, Tree}, ord::{Register, Set}, phys::{Light, Material}, sim::{Attribute, Settings}, @@ -22,8 +22,6 @@ pub struct Input<'a> { pub light: Light<'a>, /// Hit-scan tree. pub tree: &'a Tree<'a, Attribute<'a>>, - /// Measurement grid. - pub grid: &'a Grid, /// General settings. pub sett: &'a Settings, /// Boundary for the simulation. @@ -40,7 +38,6 @@ impl<'a> Input<'a> { attrs: &'a Set, light: Light<'a>, tree: &'a Tree, - grid: &'a Grid, sett: &'a Settings, bound: &'a Boundary ) -> Self { @@ -50,7 +47,6 @@ impl<'a> Input<'a> { attrs, light, tree, - grid, sett, bound, } @@ -66,7 +62,6 @@ impl Display for Input<'_> { fmt_report!(fmt, self.attrs, "attributes"); fmt_report!(fmt, self.light, "light"); fmt_report!(fmt, self.tree, "hit-scan tree"); - fmt_report!(fmt, self.grid, "measurement grid"); fmt_report!(fmt, self.sett, "settings"); fmt_report!(fmt, self.bound, "boundary"); Ok(()) diff --git a/src/sim/param/parameters.rs b/src/sim/param/parameters.rs index 07bda3a..8541d17 100644 --- a/src/sim/param/parameters.rs +++ b/src/sim/param/parameters.rs @@ -1,7 +1,7 @@ //! Runtime parameters. use crate::{ - fmt_report, geom::{Grid, SurfaceLinker, TreeSettings}, io::output::OutputConfig, ord::Set, phys::{LightLinker, Material}, sim::{AttributeLinkerLinkerLinkerLinkerLinker, Engine, Settings} + fmt_report, geom::{Boundary, SurfaceLinker, TreeSettings}, io::output::OutputConfig, ord::Set, phys::{LightLinker, Material}, sim::{AttributeLinkerLinkerLinkerLinkerLinker, Engine, Settings} }; use std::fmt::{Display, Error, Formatter}; @@ -9,10 +9,10 @@ use std::fmt::{Display, Error, Formatter}; pub struct Parameters { /// Simulation specific settings. pub sett: Settings, + /// Boundary settings. + pub boundary: Boundary, /// Tree settings. pub tree: TreeSettings, - /// Measurement grid settings. - pub grid: Grid, /// Surfaces. pub surfs: Set, /// Attributes. @@ -34,8 +34,8 @@ impl Parameters { #[inline] pub const fn new( sett: Settings, + boundary: Boundary, tree: TreeSettings, - grid: Grid, surfs: Set, attrs: Set, mats: Set, @@ -45,8 +45,8 @@ impl Parameters { ) -> Self { Self { sett, + boundary, tree, - grid, surfs, attrs, mats, @@ -62,8 +62,8 @@ impl Display for Parameters { fn fmt(&self, fmt: &mut Formatter) -> Result<(), Error> { writeln!(fmt, "...")?; fmt_report!(fmt, self.sett, "settings"); + fmt_report!(fmt, self.boundary, "boundary"); fmt_report!(fmt, self.tree, "tree settings"); - fmt_report!(fmt, self.grid, "grid settings"); fmt_report!(fmt, self.surfs, "surfaces"); fmt_report!(fmt, self.attrs, "attributes"); fmt_report!(fmt, self.mats, "materials"); diff --git a/src/sim/param/parameters_builder.rs b/src/sim/param/parameters_builder.rs index 449fa02..8d0f125 100644 --- a/src/sim/param/parameters_builder.rs +++ b/src/sim/param/parameters_builder.rs @@ -1,7 +1,7 @@ //! Buildable parameters. use crate::{ - fmt_report, geom::{GridBuilder, SurfaceLinker, TreeSettings}, io::output::OutputConfig, ord::{Build, Set}, phys::{LightLinkerBuilder, MaterialBuilder}, sim::{AttributeLinkerLinkerLinkerLinkerLinker, EngineBuilder, Parameters, Settings} + fmt_report, geom::{BoundaryBuilder, SurfaceLinker, TreeSettings}, io::output::OutputConfig, ord::{Build, Set}, phys::{LightLinkerBuilder, MaterialBuilder}, sim::{AttributeLinkerLinkerLinkerLinkerLinker, EngineBuilder, Parameters, Settings} }; use std::fmt::{Display, Error, Formatter}; @@ -9,10 +9,10 @@ use std::fmt::{Display, Error, Formatter}; pub struct ParametersBuilder { /// Simulation specific settings. sett: Settings, + /// Boundary settings. + boundary: BoundaryBuilder, /// Tree settings. tree: TreeSettings, - /// Measurement grid settings. - grid: GridBuilder, /// Surfaces. surfs: Set, /// Attributes. @@ -34,8 +34,8 @@ impl ParametersBuilder { #[inline] pub const fn new( sett: Settings, + boundary: BoundaryBuilder, tree: TreeSettings, - grid: GridBuilder, surfs: Set, attrs: Set, mats: Set, @@ -45,8 +45,8 @@ impl ParametersBuilder { ) -> Self { Self { sett, + boundary, tree, - grid, surfs, attrs, mats, @@ -63,8 +63,8 @@ impl Build for ParametersBuilder { #[inline] fn build(self) -> Self::Inst { let sett = self.sett; + let boundary = self.boundary.build(); let tree = self.tree; - let grid = self.grid.build(); let surfs = self.surfs; let attrs = self.attrs; let mats = self.mats.build(); @@ -72,7 +72,7 @@ impl Build for ParametersBuilder { let engine = self.engine.build(); let output = self.output; - Self::Inst::new(sett, tree, grid, surfs, attrs, mats, light, engine, output) + Self::Inst::new(sett, boundary, tree, surfs, attrs, mats, light, engine, output) } } @@ -81,8 +81,8 @@ impl Display for ParametersBuilder { fn fmt(&self, fmt: &mut Formatter) -> Result<(), Error> { writeln!(fmt, "...")?; fmt_report!(fmt, self.sett, "settings"); + fmt_report!(fmt, self.boundary, "boundary"); fmt_report!(fmt, self.tree, "tree settings"); - fmt_report!(fmt, self.grid, "grid settings"); fmt_report!(fmt, self.surfs, "surfaces"); fmt_report!(fmt, self.attrs, "attributes"); fmt_report!(fmt, self.mats, "materials"); diff --git a/src/sim/param/parameters_builder_loader.rs b/src/sim/param/parameters_builder_loader.rs index 596d1e6..0ca6a5d 100644 --- a/src/sim/param/parameters_builder_loader.rs +++ b/src/sim/param/parameters_builder_loader.rs @@ -1,7 +1,7 @@ //! Loadable parameters. use crate::{ - err::Error, fs::{Load, Redirect}, geom::{GridBuilder, SurfaceLinkerLoader, TreeSettings}, io::output::OutputConfig, ord::Set, phys::{LightLinkerBuilderLoader, MaterialBuilder}, sim::{ + err::Error, fs::{Load, Redirect}, geom::{boundary_builder::BoundaryBuilder, SurfaceLinkerLoader, TreeSettings}, io::output::OutputConfig, ord::Set, phys::{LightLinkerBuilderLoader, MaterialBuilder}, sim::{ AttributeLinkerLinkerLinkerLinkerLinker, EngineBuilderLoader, ParametersBuilder, Settings, } }; @@ -13,10 +13,10 @@ use std::path::Path; pub struct ParametersBuilderLoader { /// Simulation specific settings. sett: Redirect, + // Boundary conditions. + boundary: Redirect, /// Tree settings. tree: Redirect, - /// Measurement grid settings. - grid: Redirect, /// Surfaces. surfs: Redirect>, /// Attributes. @@ -37,8 +37,8 @@ impl Load for ParametersBuilderLoader { #[inline] fn load(self, in_dir: &Path) -> Result { let sett = self.sett.load(in_dir)?; + let boundary = self.boundary.load(in_dir)?; let tree = self.tree.load(in_dir)?; - let grid = self.grid.load(in_dir)?; let surfs = self.surfs.load(in_dir)?.load(in_dir)?; let attrs = self.attrs.load(in_dir)?; let mats = self.mats.load(in_dir)?.load(in_dir)?; @@ -47,7 +47,7 @@ impl Load for ParametersBuilderLoader { let output = self.output.load(in_dir)?; Ok(Self::Inst::new( - sett, tree, grid, surfs, attrs, mats, lights, engine, output, + sett, boundary, tree, surfs, attrs, mats, lights, engine, output, )) } } From bf69dc2636a410c2b629ebffd7b542b20ed28e73 Mon Sep 17 00:00:00 2001 From: Sam Morrell Date: Thu, 15 Aug 2024 23:25:34 +0100 Subject: [PATCH 16/16] Added potential fix for shift map indexing in `fluorescence` engine. The removal of the `grid` parameter also removed access to a consistent measurement grid over which to index into the shifts map in the fluorescence engine. However, I have engineered an alternate solution which indexes using the bounds of the simulation and takes from the res of the grid, potentially making this solution more flexible. --- src/geom/domain/boundary.rs | 23 ++++++++++++++++++++++- src/sim/engine/engines/fluorescence.rs | 6 +++--- 2 files changed, 25 insertions(+), 4 deletions(-) diff --git a/src/geom/domain/boundary.rs b/src/geom/domain/boundary.rs index 024ec5f..7c8d985 100644 --- a/src/geom/domain/boundary.rs +++ b/src/geom/domain/boundary.rs @@ -3,7 +3,8 @@ use crate::{ geom::{plane::ray_plane_intersection, Cube, Hit, Ray, Side, Trace}, math::{Dir3, Point3, Vec3}, phys::{Photon, Reflectance}, - sim::Attribute + sim::Attribute, + ord::cartesian::{X, Y, Z}, }; use rand::rngs::ThreadRng; use std::fmt::{Display, Formatter}; @@ -229,6 +230,26 @@ impl Boundary { return facing_dir; } + + /// If the given position is contained within the grid, + /// generate the index for the given position within the grid. + #[inline] + #[must_use] + pub fn gen_index(&self, p: &Point3, res: [usize; 3]) -> Option<[usize; 3]> { + self.contains(p).then(|| { + let mins = self.bounding_box.mins(); + let maxs = self.bounding_box.maxs(); + + [ + (((p.x() - mins.x()) / (maxs.x() - mins.x())) * res[X] as f64).floor() + as usize, + (((p.y() - mins.y()) / (maxs.y() - mins.y())) * res[Y] as f64).floor() + as usize, + (((p.z() - mins.z()) / (maxs.z() - mins.z())) * res[Z] as f64).floor() + as usize, + ] + }) + } } impl Display for Boundary { diff --git a/src/sim/engine/engines/fluorescence.rs b/src/sim/engine/engines/fluorescence.rs index 5864f7d..29f0899 100644 --- a/src/sim/engine/engines/fluorescence.rs +++ b/src/sim/engine/engines/fluorescence.rs @@ -39,6 +39,7 @@ pub fn fluorescence( let mat = input.light.mat(); let mut local = mat.sample_environment(phot.wavelength()); let mut env; + let flu_concs_res = [flu_concs.shape()[0], flu_concs.shape()[1], flu_concs.shape()[2]]; // Main event loop. let mut num_loops = 0; @@ -60,9 +61,8 @@ pub fn fluorescence( } // Local variable modifications. - // TODO: I have had to remove this for now, as I've removed the fixed grid. - //let index = input.grid.gen_index_voxel(phot.ray().pos()); - let index = [0, 0, 0]; + let index = input.bound.gen_index(phot.ray().pos(), flu_concs_res) + .expect("Unable to index shift map with photon position. "); env = Local::new( local.ref_index(), local.scat_coeff(),