diff --git a/src/mat.rs b/src/mat.rs index 6c4bcb3..e42fc62 100644 --- a/src/mat.rs +++ b/src/mat.rs @@ -1,5 +1,4 @@ use crate::CellState; -use crate::flt::{self, Flt}; use crate::geom::{Line2d, Vec2, Vec3, Polygon2d}; use crate::real::Real; use crate::sim::{PmlParameters, PmlState, StepParameters, StepParametersMut}; diff --git a/src/post.rs b/src/post.rs index 4214807..9ce505f 100644 --- a/src/post.rs +++ b/src/post.rs @@ -1,5 +1,5 @@ //! Post-processing tools -use crate::mat::GenericMaterial; +use crate::mat::{GenericMaterial, Static}; use crate::meas::AbstractMeasurement; use crate::render::{ColorTermRenderer, Renderer as _, RenderConfig, SerializedFrame}; use crate::sim::{SimState, StaticSim}; @@ -34,7 +34,7 @@ pub type Result = std::result::Result; pub struct Frame { path: PathBuf, - data: SerializedFrame, + data: SerializedFrame, } impl Frame { @@ -111,9 +111,9 @@ impl Loader { // decode to a valid but incorrect state... let data = bincode::deserialize_from(&mut reader).or_else(|_| -> Result<_> { reader.seek(SeekFrom::Start(0)).unwrap(); - let data: SerializedFrame>> = + let data: SerializedFrame, f32>> = bincode::deserialize_from(reader)?; - Ok(data.to_static()) + Ok(SerializedFrame::to_static(data)) })?; Ok(Frame { path: path.into(), diff --git a/src/real.rs b/src/real.rs index 59c2767..e9e67ce 100644 --- a/src/real.rs +++ b/src/real.rs @@ -4,6 +4,16 @@ use decorum::cmp::IntrinsicOrd; pub use decorum::Real as decorum_Real; pub use num::Zero; +// finite. non-nan, non-inf +pub mod checked { + pub use decorum::{R32, R64}; +} + +pub mod unchecked { + pub type R32 = f32; + pub type R64 = f64; +} + pub trait ToFloat { fn to_f32(&self) -> f32 { self.to_f64() as _ diff --git a/src/render.rs b/src/render.rs index 9f2573c..94c4264 100644 --- a/src/render.rs +++ b/src/render.rs @@ -643,7 +643,7 @@ pub struct SerializedFrame { impl SerializedFrame { pub fn to_static(self) -> SerializedFrame { SerializedFrame { - state: self.state.to_static(), + state: GenericSim::to_static(&self.state), measurements: self.measurements, } } diff --git a/src/sim.rs b/src/sim.rs index ee3e9d7..761f5d6 100644 --- a/src/sim.rs +++ b/src/sim.rs @@ -1,9 +1,7 @@ -use crate::flt::Real; use crate::geom::{Coord, Cube, Index, InvertedRegion, Meters, Region, Vec3, Vec3u}; use crate::mat::{self, GenericMaterial, Material, MaterialExt as _}; -use crate::real::{self, Real as _, ToFloat as _}; +use crate::real::{checked, Real, ToFloat as _}; use crate::stim::AbstractStimulus; -use dyn_clone::{self, DynClone}; use log::trace; use ndarray::{Array3, Zip}; use rayon::prelude::*; @@ -11,7 +9,7 @@ use serde::{Serialize, Deserialize}; use std::convert::From; use std::iter::Sum; -pub type StaticSim = SimState, Real>; +pub type StaticSim = SimState, f32>; #[derive(Default, Copy, Clone, PartialEq, Serialize, Deserialize)] pub struct PmlState { @@ -46,7 +44,7 @@ struct PmlAux { #[derive(Default, Copy, Clone, PartialEq, Serialize, Deserialize)] struct ConvFlt(R); -impl ConvFlt { +impl ConvFlt { ///```tex /// Consider $v(t) = (a \delta(t) + b exp[-c t] u(t)) \ast y(t)$ /// This method will advance from self = $v(t)$ to self = $v(t + \Delta t)$ @@ -114,8 +112,8 @@ pub struct StepParametersMut<'a, R> { pml: Option<(&'a mut PmlState, PmlParameters)>, } -impl<'a, R: real::Real> StepParametersMut<'a, R> { - pub fn new( +impl<'a, R: Real> StepParametersMut<'a, R> { + pub fn new( conductivity: Vec3, pml: Option<(&'a mut PmlState, PmlParameters)> ) -> Self { Self { @@ -123,7 +121,7 @@ impl<'a, R: real::Real> StepParametersMut<'a, R> { pml, } } - pub fn with_conductivity(mut self, conductivity: Vec3) -> Self { + pub fn with_conductivity(mut self, conductivity: Vec3) -> Self { self.conductivity = conductivity.cast(); self } @@ -167,8 +165,8 @@ pub struct PmlParameters { pseudo_conductivity: Vec3, } -impl PmlParameters { - pub fn new(pseudo_conductivity: Vec3) -> Self { +impl PmlParameters { + pub fn new(pseudo_conductivity: Vec3) -> Self { Self { pseudo_conductivity: pseudo_conductivity.cast(), } @@ -182,7 +180,7 @@ pub struct StepParameters<'a, R> { pml: Option<(&'a PmlState, PmlParameters)>, } -impl<'a, R: real::Real> StepParameters<'a, R> { +impl<'a, R: Real> StepParameters<'a, R> { pub fn conductivity(&self) -> Vec3 { self.conductivity } @@ -200,7 +198,7 @@ impl<'a, R> From> for StepParameters<'a, R> { } } -pub trait GenericSim: Send + Sync + DynClone { +pub trait GenericSim: Send + Sync { fn sample(&self, pos: Meters) -> Sample; /// DEPRECATED. Use stimulus instead fn impulse_e_meters(&mut self, pos: Meters, amount: Vec3); @@ -257,7 +255,6 @@ pub trait GenericSim: Send + Sync + DynClone { state } } -dyn_clone::clone_trait_object!(GenericSim); impl<'a> dyn GenericSim + 'a { pub fn get(&self, at: C) -> Sample { @@ -372,7 +369,7 @@ impl<'a> dyn GenericSim + 'a { } #[derive(Default, Clone, Serialize, Deserialize)] -pub struct SimState, R=Real> { +pub struct SimState, R=checked::R32> { cells: Array3, e: Array3>, h: Array3>, @@ -400,7 +397,7 @@ impl SimState { } } -impl SimState { +impl SimState { pub fn feature_size(&self) -> f32 { self.feature_size.to_f32() } @@ -409,7 +406,7 @@ impl SimState { } } -impl SimState { +impl SimState { pub fn new(size: Index, feature_size: f32) -> Self { // XXX this has to be multiplied by a Courant Number in order to ensure stability. // For 3d, we need 3 full timesteps in order to communicate information across the corner @@ -431,7 +428,7 @@ impl SimState { } } -impl + Clone + Send + Sync + 'static> SimState { +impl + Send + Sync> SimState { pub fn step(&mut self) { let time_step = self.timestep; let half_time_step = R::half() * time_step; @@ -479,13 +476,13 @@ impl + Clone + Send + Sync + 'static> SimState + Clone> SimState { +impl + Clone> SimState { pub fn fill_region(&mut self, region: &Reg, mat: M) { self.fill_region_using(region, |_idx: Index| mat.clone()); } } -impl> SimState { +impl> SimState { pub fn fill_region_using(&mut self, region: &Reg, f: F) where Reg: Region, @@ -554,7 +551,7 @@ impl> SimState { } } -impl + Clone + Send + Sync + 'static> GenericSim for SimState { +impl + Send + Sync> GenericSim for SimState { fn sample(&self, pos: Meters) -> Sample { // TODO: smarter sampling than nearest neighbor? let pos_sim = pos.to_index(self.feature_size()); @@ -610,7 +607,7 @@ impl + Clone + Send + Sync + 'static> GenericSim f } #[allow(unused)] -impl SimState { +impl SimState { fn get_mat(&self, c: C) -> &M { let at = c.to_index(self.feature_size.cast()); &self.cells[at.row_major_idx()] @@ -800,7 +797,7 @@ impl Sample { /// Solve df/dt in: /// nabla_g = df_dt_coeff * df/dt + f_coeff * f + f_int_coeff * \int f + f_int_int_coeff * \int \int f #[allow(unused)] -fn solve_step_diff_eq( +fn solve_step_diff_eq( nabla_g: Vec3, df_dt_coeff: Vec3, f_coeff: Vec3, @@ -840,7 +837,7 @@ fn solve_step_diff_eq( /// Solve df/dt in: /// nabla_g = df_dt_coeff * df/dt + f_coeff * f /// Lower-order (more efficient) version of `solve_step_diff_eq`. -fn solve_step_diff_eq_linear( +fn solve_step_diff_eq_linear( nabla_g: Vec3, df_dt_coeff: Vec3, f_coeff: Vec3, @@ -868,7 +865,7 @@ struct ECell<'a, M, R> { h: &'a Vec3, } -impl<'a, R: real::Real, M: Material> HCell<'a, M, R> { +impl<'a, R: Real, M: Material> HCell<'a, M, R> { pub fn b(&self) -> Vec3 { (*self.h + self.mat.m()) * R::mu0() } @@ -1039,7 +1036,7 @@ impl<'a, R: real::Real, M: Material> HCell<'a, M, R> { } } -impl<'a, R: real::Real, M: Material> ECell<'a, M, R> { +impl<'a, R: Real, M: Material> ECell<'a, M, R> { /// delta_t = timestep covered by this function. i.e. it should be half the timestep of the simulation /// feature_size = how many units apart is the center of each adjacent cell on the grid. fn step(&mut self, neighbors: Neighbors<'_, R>, delta_t: R, inv_feature_size: R) { @@ -1146,7 +1143,7 @@ pub struct CellState { h: Vec3, } -impl CellState { +impl CellState { pub fn h(&self) -> Vec3 { self.h } @@ -1281,7 +1278,7 @@ mod test { do_conv_flt_test(a, b, c, delta_t, &signal, &expected); } - fn assert_vec3_eq(actual: Vec3, expected: Vec3) { + fn assert_vec3_eq(actual: Vec3, expected: Vec3) { let actual = actual.cast::(); let expected = expected.cast::(); assert_float_eq!(actual.x(), expected.x(), r2nd <= 1e-6); @@ -1521,7 +1518,7 @@ mod test { #[test] fn accessors() { - let mut state = StaticSim::new(Index((15, 17, 19).into()), 1e-6); + let mut state = SimState::, checked::R64>::new(Index((15, 17, 19).into()), 1e-6); assert_eq!(state.width(), 15); assert_eq!(state.height(), 17); assert_eq!(state.depth(), 19); @@ -1534,7 +1531,7 @@ mod test { #[test] fn energy_conservation_over_time() { - let mut state = StaticSim::new(Index((2001, 1, 1).into()), 1e-6); + let mut state = SimState::, checked::R64>::new(Index((2001, 1, 1).into()), 1e-6); for t in 0..100 { let angle = (t as f32)/100.0*f32::two_pi(); let sig = angle.sin(); @@ -1552,7 +1549,7 @@ mod test { #[test] fn sane_boundary_conditions() { - let mut state = StaticSim::new(Index((21, 21, 21).into()), 1e-6); + let mut state = SimState::, checked::R64>::new(Index((21, 21, 21).into()), 1e-6); for t in 0..10 { let angle = (t as f32)/10.0*f32::two_pi(); let sig = angle.sin(); @@ -1569,8 +1566,8 @@ mod test { /// Fill the world with the provided material and a stimulus. /// Measure energy at the start, and then again after advancing many steps. /// Return these two measurements (energy(t=0), energy(t=~=1000)) - fn conductor_test>>(mat: M) -> (f32, f32) { - let mut state = StaticSim::new(Index((201, 1, 1).into()), 1e-6); + fn conductor_test>>(mat: M) -> (f32, f32) { + let mut state = SimState::, checked::R64>::new(Index((201, 1, 1).into()), 1e-6); state.fill_region(&WorldRegion, mat.into()); for t in 0..100 { let angle = (t as f32)/100.0*f32::two_pi(); @@ -1616,7 +1613,7 @@ mod test { assert_float_eq!(energy_1/energy_0, 0.0, abs <= 1e-6); } - fn state_for_pml(size: Index) -> SimState> { + fn state_for_pml(size: Index) -> SimState, checked::R64> { let mut state = SimState::new(size, 1e-6); let timestep = state.timestep(); state.fill_boundary_using(size/4, |boundary_ness| { @@ -1635,15 +1632,13 @@ mod test { } /// Like state_for_pml, but the boundary is an ordinary electric conductor -- not PML - fn state_for_graded_conductor(size: Index) -> StaticSim { - let mut state = StaticSim::new(size, 1e-6); + fn state_for_graded_conductor(size: Index) -> SimState, checked::R64> { + let mut state = SimState::new(size, 1e-6); let timestep = state.timestep(); state.fill_boundary_using(size/4, |boundary_ness| { let b = boundary_ness.elem_pow(3.0); let conductivity = Vec3::uniform(f32::eps0() * b.mag() * 0.5 / timestep); - Static { - conductivity: conductivity.cast(), ..Default::default() - } + Conductor::new_anisotropic(conductivity) }); state } @@ -1782,7 +1777,7 @@ mod test { pml_test_against_baseline(&mut pml_state, &mut baseline_state, Vec3::unit_x()); } - fn state_for_monodirectional_pml(size: Index) -> SimState> { + fn state_for_monodirectional_pml(size: Index) -> SimState, checked::R64> { let mut state = SimState::new(size, 1e-6); let timestep = state.timestep(); state.fill_boundary_using(size/4, |boundary_ness| { @@ -1805,7 +1800,7 @@ mod test { fn pml_ineffective_mono_linear_test) -> Vec3>( size: Index, e: Vec3, shuffle: F ) { - let mut state = SimState::>::new(size, 1e-6); + let mut state = SimState::new(size, 1e-6); let timestep = state.timestep(); state.fill_boundary_using(size/4, |boundary_ness| { let b = boundary_ness.elem_pow(3.0); @@ -1813,7 +1808,7 @@ mod test { // let cs = b * 0.5; // permute the stretching so that it shouldn't interfere with the wave let coord_stretch = shuffle(cs); - Pml::new(coord_stretch) + Pml::::new(coord_stretch) }); let center = Index(state.size().0 / 2); let en = pml_test_at(&mut state, e, center);