From 8779371f427d39e33d134cbab8e935c8b8db68c7 Mon Sep 17 00:00:00 2001 From: Colin Date: Mon, 7 Jun 2021 18:25:11 -0700 Subject: [PATCH] Stabilize time as f32. --- examples/minimal_torus.rs | 16 ++--- examples/toroid25d.rs | 12 ++-- examples/wrapped_torus.rs | 22 +++--- src/bin/pml_tuning.rs | 39 ++++++----- src/driver.rs | 18 ++--- src/geom/region.rs | 56 +++++++-------- src/geom/units.rs | 37 +++++----- src/geom/vecu.rs | 6 +- src/meas.rs | 75 ++++++++++---------- src/real.rs | 58 ++++++++++++---- src/sim.rs | 143 +++++++++++++++++++------------------- src/stim.rs | 67 +++++++++--------- 12 files changed, 288 insertions(+), 261 deletions(-) diff --git a/examples/minimal_torus.rs b/examples/minimal_torus.rs index 58e1f88..fa01973 100644 --- a/examples/minimal_torus.rs +++ b/examples/minimal_torus.rs @@ -4,7 +4,7 @@ use coremem::stim::{CurlStimulus, Sinusoid1, TimeVarying1 as _}; fn main() { coremem::init_logging(); - let feat_size = 10e-6; // feature size + let feat_size = 10e-6f32; // feature size let duration = 6e-9; let width = 2200e-6; let depth = 1800e-6; @@ -16,10 +16,10 @@ fn main() { let peak_current = 5e11; let current_duration = 1.0e-11; // half-wavelength of the sine wave let current_break = 0.5e-11; // time between 'set' pulse and 'clear' pulse - let conductivity = 1.0e9; + let conductivity = 1.0e9f32; - let from_m = |m: Flt| (m/feat_size).round() as u32; - let m_to_um = |m: Flt| (m * 1e6).round() as u32; + let from_m = |m: f32| (m/feat_size).round() as u32; + let m_to_um = |m: f32| (m * 1e6).round() as u32; let half_width = width * 0.5; let half_depth = depth * 0.5; @@ -49,8 +49,8 @@ fn main() { ); driver.fill_region(&ferro_region, mat::db::minimal_square_ferrite().into()); - driver.fill_region(&drive_region, mat::db::conductor(conductivity).into()); - driver.fill_region(&sense_region, mat::db::conductor(conductivity).into()); + driver.fill_region(&drive_region, mat::db::conductor(conductivity as _).into()); + driver.fill_region(&sense_region, mat::db::conductor(conductivity as _).into()); // driver.add_pml_boundary(Meters((boundary_xy, boundary_xy, boundary_z).into())); driver.add_classical_boundary(Meters::new(boundary_xy, boundary_xy, boundary_z)); @@ -61,10 +61,10 @@ fn main() { // if I = k*sin(w t) then dE/dt = k*w sin(w t) / (A*\sigma) // i.e. dE/dt is proportional to I/(A*\sigma), multiplied by w (or, divided by wavelength) let peak_stim = peak_current/current_duration / (drive_region.cross_section() * conductivity); - let pos_wave = Sinusoid1::from_wavelength(peak_stim, current_duration * 2.0) + let pos_wave = Sinusoid1::from_wavelength(peak_stim as _, current_duration * 2.0) .half_cycle(); - let neg_wave = Sinusoid1::from_wavelength(-peak_stim, current_duration * 2.0) + let neg_wave = Sinusoid1::from_wavelength(-peak_stim as _, current_duration * 2.0) .half_cycle() .shifted(current_duration + current_break); diff --git a/examples/toroid25d.rs b/examples/toroid25d.rs index 34ffe1f..ce2951f 100644 --- a/examples/toroid25d.rs +++ b/examples/toroid25d.rs @@ -6,11 +6,11 @@ use log::trace; fn main() { coremem::init_logging(); for p in 0..20 { - let ferro_depth = 10e-6 * (20-p) as Flt; - let feat_size = 10e-6; // feature size + let ferro_depth = 10e-6 * (20-p) as f32; + let feat_size = 10e-6f32; // feature size let from_m = |m| (m/feat_size) as u32; let m_to_um = |px| (px * 1e6) as u32; - let to_m = |px| px as Flt * feat_size; + let to_m = |px| px as f32 * feat_size; let width = 2500e-6; let depth = 200e-6; let conductor_inner_rad = 0e-6; @@ -91,7 +91,7 @@ fn main() { Meters::new(half_width + 0.5 * (ferro_inner_rad + ferro_outer_rad), half_width, half_depth) )); - let center = Vec2::new(half_width as Flt, half_width as Flt); + let center = Vec2::new(half_width, half_width); for z_px in 0..depth_px { for y_px in 0..width_px { @@ -99,9 +99,9 @@ fn main() { let loc = Index((x_px, y_px, z_px).into()); let d = Vec2::new(to_m(x_px), to_m(y_px)) - center; let r = d.mag(); - if (conductor_inner_rad as Flt..conductor_outer_rad as Flt).contains(&r) { + if (conductor_inner_rad..conductor_outer_rad).contains(&r) { driver.state.put_material(loc, mat::db::conductor(conductivity).into()); - } else if (ferro_inner_rad as Flt..ferro_outer_rad as Flt).contains(&r) { + } else if (ferro_inner_rad..ferro_outer_rad).contains(&r) { let half_depth_px = from_m(half_depth); let ferro_depth_px = from_m(ferro_depth); if (half_depth_px-ferro_depth_px/2 .. half_depth_px+(ferro_depth_px+1)/2).contains(&z_px) { diff --git a/examples/wrapped_torus.rs b/examples/wrapped_torus.rs index dfe562d..4f806cc 100644 --- a/examples/wrapped_torus.rs +++ b/examples/wrapped_torus.rs @@ -4,7 +4,7 @@ use coremem::stim::{CurlStimulus, Sinusoid1, TimeVarying1 as _}; fn main() { coremem::init_logging(); - let feat_size = 10e-6; // feature size + let feat_size = 10e-6f32; // feature size let duration = 0.5e-9; let width = 4800e-6; let height = 2200e-6; @@ -17,11 +17,11 @@ fn main() { let peak_current = 5e6; let current_duration = 0.1e-9; // half-wavelength of the sine wave let current_break = 0.1e-9; // time between 'set' pulse and 'clear' pulse - let drive_conductivity = 5.0; - let sense_conductivity = 5.0; + let drive_conductivity = 5.0f32; + let sense_conductivity = 5.0f32; - let from_m = |m: Flt| (m/feat_size).round() as u32; - let m_to_um = |m: Flt| (m * 1e6).round() as u32; + let from_m = |m: f32| (m/feat_size).round() as u32; + let m_to_um = |m: f32| (m * 1e6).round() as u32; let half_width = width * 0.5; let q1_width = width * 0.25; let q3_width = width * 0.75; @@ -42,16 +42,16 @@ fn main() { let sense1_region = Torus::new_xz(Meters::new(q1_width + ferro_major, half_height, half_depth), wire_major, wire_minor); driver.fill_region(&ferro1_region, mat::db::linear_iron().into()); - driver.fill_region(&drive1_region, mat::db::conductor(drive_conductivity).into()); - driver.fill_region(&sense1_region, mat::db::conductor(sense_conductivity).into()); + driver.fill_region(&drive1_region, mat::db::conductor(drive_conductivity as _).into()); + driver.fill_region(&sense1_region, mat::db::conductor(sense_conductivity as _).into()); let ferro2_region = Torus::new_xy(Meters::new(q3_width, half_height, half_depth), ferro_major, ferro_minor); let drive2_region = Torus::new_xz(Meters::new(q3_width - ferro_major, half_height, half_depth), wire_major, wire_minor); let sense2_region = Torus::new_xz(Meters::new(q3_width + ferro_major, half_height, half_depth), wire_major, wire_minor); driver.fill_region(&ferro2_region, mat::db::ferroxcube_3r1().into()); - driver.fill_region(&drive2_region, mat::db::conductor(drive_conductivity).into()); - driver.fill_region(&sense2_region, mat::db::conductor(sense_conductivity).into()); + driver.fill_region(&drive2_region, mat::db::conductor(drive_conductivity as _).into()); + driver.fill_region(&sense2_region, mat::db::conductor(sense_conductivity as _).into()); let boundary_xy = q1_width - ferro_major - ferro_minor - buffer; let boundary_z = half_depth - wire_major - wire_minor - buffer; @@ -65,10 +65,10 @@ fn main() { // if I = k*sin(w t) then dE/dt = k*w sin(w t) / (A*\sigma) // i.e. dE/dt is proportional to I/(A*\sigma), multiplied by w (or, divided by wavelength) let peak_stim = peak_current/current_duration / (drive1_region.cross_section() * drive_conductivity); - let pos_wave = Sinusoid1::from_wavelength(peak_stim, current_duration * 2.0) + let pos_wave = Sinusoid1::from_wavelength(peak_stim as Flt, current_duration * 2.0) .half_cycle(); - let neg_wave = Sinusoid1::from_wavelength(-4.0*peak_stim, current_duration * 2.0) + let neg_wave = Sinusoid1::from_wavelength(-4.0*peak_stim as Flt, current_duration * 2.0) .half_cycle() .shifted(current_duration + current_break); diff --git a/src/bin/pml_tuning.rs b/src/bin/pml_tuning.rs index 8bad120..bd02d6d 100644 --- a/src/bin/pml_tuning.rs +++ b/src/bin/pml_tuning.rs @@ -1,18 +1,18 @@ use coremem::*; use coremem::geom::*; -use coremem::real::Real as _; +use coremem::real::{Real as _, ToFloat as _}; use coremem::stim::AbstractStimulus; use rand::rngs::StdRng; use rand::{Rng as _, SeedableRng as _}; -fn energy(s: &dyn GenericSim, reg: &R) -> Flt { - let e = Real::half() * s.map_sum_over_enumerated(reg, |_pos: Index, cell| { - cell.e().mag_sq() +fn energy(s: &dyn GenericSim, reg: &R) -> f32 { + let e = f64::half() * s.map_sum_over_enumerated(reg, |_pos: Index, cell| { + cell.e().mag_sq().to_f64() }); - Flt::from_primitive(e) + e.cast() } -fn energy_now_and_then(state: &mut StaticSim, reg: &R, frames: u32) -> (Flt, Flt) { +fn energy_now_and_then(state: &mut StaticSim, reg: &R, frames: u32) -> (f32, f32) { let energy_0 = energy(state, reg); for _ in 0..frames { state.step(); @@ -22,13 +22,14 @@ fn energy_now_and_then(state: &mut StaticSim, reg: &R, frames: u32) - } struct PmlStim { + /// Maps index -> (stim vector, stim frequency) f: F, - t_end: Flt, - feat_size: Flt, + t_end: f32, + feat_size: f32, } -impl (Vec3, Flt) + Sync> AbstractStimulus for PmlStim { - fn at(&self, t_sec: Flt, pos: Meters) -> Vec3 { - let angle = (t_sec as Flt)/self.t_end*2.0*consts::PI; +impl (Vec3, f32) + Sync> AbstractStimulus for PmlStim { + fn at(&self, t_sec: f32, pos: Meters) -> Vec3 { + let angle = t_sec/self.t_end*f32::two_pi(); let gate = 0.5*(1.0 - angle.cos()); let (e, hz) = (self.f)(pos.to_index(self.feat_size)); let sig_angle = angle*hz; @@ -39,11 +40,11 @@ impl (Vec3, Flt) + Sync> AbstractStimulus for PmlStim { /// Apply some stimulus, and then let it decay and measure the ratio of energy left in the system fn apply_stim_full_interior(state: &mut StaticSim, frames: u32, f: F) -where F: Fn(Index) -> (Vec3, Flt) + Sync // returns (E vector, omega) +where F: Fn(Index) -> (Vec3, f32) + Sync // returns (E vector, omega) { let stim = PmlStim { f, - t_end: (frames as Flt) * state.timestep(), + t_end: (frames as f32) * state.timestep(), feat_size: state.feature_size(), }; for _t in 0..frames { @@ -55,7 +56,7 @@ where F: Fn(Index) -> (Vec3, Flt) + Sync // returns (E vector, omega) fn apply_stim_over_region(state: &mut StaticSim, frames: u32, reg: R, f: F) where R: Region, - F: Fn(Index) -> (Vec3, Flt) + Sync, + F: Fn(Index) -> (Vec3, f32) + Sync, { let feat = state.feature_size(); apply_stim_full_interior(state, frames, |idx| { @@ -84,7 +85,7 @@ fn apply_chaotic_stim_over_region(state: &mut StaticSim, frames: u32, }) } -fn chaotic_pml_test(state: &mut StaticSim, boundary: u32, padding: u32, frames: u32) -> Flt { +fn chaotic_pml_test(state: &mut StaticSim, boundary: u32, padding: u32, frames: u32) -> f32 { let feat = state.feature_size(); { let upper_left_idx = Index::unit()*padding; @@ -102,12 +103,12 @@ fn chaotic_pml_test(state: &mut StaticSim, boundary: u32, padding: u32, frames: energy_1/energy_0 } -fn state_for_pml(size: Index, boundary: Index, feat_size: Flt, sc_coeff: Flt, cond_coeff: Flt, pow: Flt) -> StaticSim { +fn state_for_pml(size: Index, boundary: Index, feat_size: f32, sc_coeff: f32, cond_coeff: f32, pow: f32) -> StaticSim { let mut state = StaticSim::new(size, feat_size); let timestep = state.timestep(); state.fill_boundary_using(boundary, |boundary_ness| { - let b = boundary_ness.elem_pow(Real::from_primitive(pow)); - let coord_stretch = b * Real::from_primitive(sc_coeff / timestep); + let b = boundary_ness.elem_pow(pow); + let coord_stretch = b * sc_coeff / timestep; let conductivity = Vec3::unit() * (b.mag() * cond_coeff / timestep); unimplemented!() // Static { @@ -168,7 +169,7 @@ fn main() { 0.5, ] { //for cond_coeff in vec![0.0, 1.0, 1e3, 1e6, 1e9] { - for cond_coeff in vec![0.0, 0.5*consts::EPS0] { + for cond_coeff in vec![0.0, 0.5*f32::eps0()] { for frames in vec![400, 1200] { for pml_width in vec![1, 2, 4, 8, 16] { for feat_size in vec![1e-6] { diff --git a/src/driver.rs b/src/driver.rs index 3c05f66..20ef19b 100644 --- a/src/driver.rs +++ b/src/driver.rs @@ -34,7 +34,7 @@ pub struct Driver { } impl Driver { - pub fn new(size: C, feature_size: Flt) -> Self { + pub fn new(size: C, feature_size: f32) -> Self { Driver { state: SimState::new(size.to_index(feature_size), feature_size), renderer: Arc::new(MultiRenderer::new()), @@ -180,7 +180,7 @@ impl Driver { } } - pub fn step_until(&mut self, deadline: Flt) { + pub fn step_until(&mut self, deadline: f32) { while self.dyn_state().time() < deadline { self.step(); } @@ -194,9 +194,9 @@ impl> Driver { pub fn add_pml_boundary(&mut self, thickness: C) { let timestep = self.state.timestep(); self.state.fill_boundary_using(thickness, |boundary_ness| { - let b = boundary_ness.elem_pow(flt::Real::three()); - let conductivity = b * (flt::Real::half() / timestep); - Pml::new(conductivity).into() + let b = boundary_ness.elem_pow(3.0); + let conductivity = b * (0.5 / timestep); + Pml::new(conductivity.cast()).into() }); } } @@ -205,8 +205,8 @@ impl> Driver { pub fn add_classical_boundary(&mut self, thickness: C) { let timestep = self.state.timestep(); self.state.fill_boundary_using(thickness, |boundary_ness| { - let b = boundary_ness.elem_pow(flt::Real::three()); - let cond = b * (flt::Real::half() / timestep); + let b = boundary_ness.elem_pow(3.0); + let cond = b * (0.5 / timestep); let iso_cond = cond.x() + cond.y() + cond.z(); let iso_conductor = mat::db::conductor(Flt::from_primitive(iso_cond)); @@ -223,8 +223,8 @@ struct StimuliAdapter { } impl AbstractStimulus for StimuliAdapter { - fn at(&self, t_sec: Flt, pos: Meters) -> Vec3 { - self.stim.at(t_sec, pos) * flt::Real::from_f64(self.frame_interval as f64) + fn at(&self, t_sec: f32, pos: Meters) -> Vec3 { + self.stim.at(t_sec, pos) * flt::Real::from_primitive(self.frame_interval) } } diff --git a/src/geom/region.rs b/src/geom/region.rs index 2886381..752ad2d 100644 --- a/src/geom/region.rs +++ b/src/geom/region.rs @@ -1,4 +1,3 @@ -use crate::flt::{Flt, Real}; use crate::geom::{Meters, Vec2, Vec3}; use crate::real::Real as _; @@ -16,16 +15,16 @@ dyn_clone::clone_trait_object!(Region); #[derive(Copy, Clone, Serialize, Deserialize)] pub struct CylinderZ { - center: Vec2, - radius: Real, + center: Vec2, + radius: f32, } impl CylinderZ { // TODO: should be unit-typed - pub fn new(center: Vec2, radius: Flt) -> Self { + pub fn new(center: Vec2, radius: f32) -> Self { Self { center, - radius: radius.into() + radius, } } } @@ -49,24 +48,24 @@ pub struct Torus { /// Unit-length vector describing the axis of the torus normal: Meters, /// Distance from origin to the "center" of the solid part of the torus - major_rad: Real, + major_rad: f32, /// Distance from a center-point of the torus to its surface - minor_rad: Real, + minor_rad: f32, } impl Torus { - pub fn new(center: Meters, normal: Meters, major_rad: Flt, minor_rad: Flt) -> Self { + pub fn new(center: Meters, normal: Meters, major_rad: f32, minor_rad: f32) -> Self { Self { center, normal, - major_rad: Real::from_inner(major_rad), - minor_rad: Real::from_inner(minor_rad), + major_rad, + minor_rad, } } - pub fn new_xy(center: Meters, major_rad: Flt, minor_rad: Flt) -> Self { + pub fn new_xy(center: Meters, major_rad: f32, minor_rad: f32) -> Self { Self::new(center, Meters(Vec3::new(0.0, 0.0, 1.0)), major_rad, minor_rad) } - pub fn new_xz(center: Meters, major_rad: Flt, minor_rad: Flt) -> Self { + pub fn new_xz(center: Meters, major_rad: f32, minor_rad: f32) -> Self { Self::new(center, Meters(Vec3::new(0.0, 1.0, 0.0)), major_rad, minor_rad) } pub fn center(&self) -> Meters { @@ -75,9 +74,8 @@ impl Torus { pub fn axis(&self) -> Meters { self.normal } - pub fn cross_section(&self) -> Flt { - use crate::consts::PI; - (self.minor_rad * self.minor_rad * PI).into() + pub fn cross_section(&self) -> f32 { + self.minor_rad * self.minor_rad * f32::pi() } } @@ -102,21 +100,21 @@ impl Region for Torus { p_on_plane.with_mag(self.major_rad) }; let distance_to_circle = (rel_p - q).mag(); - distance_to_circle < self.minor_rad.into_inner() + distance_to_circle < self.minor_rad } } #[derive(Copy, Clone, Serialize, Deserialize)] pub struct Sphere { center: Meters, - rad: Real, + rad: f32, } impl Sphere { - pub fn new(center: Meters, rad: Flt) -> Self { + pub fn new(center: Meters, rad: f32) -> Self { Self { center, - rad: Real::from_inner(rad), + rad, } } } @@ -124,7 +122,7 @@ impl Sphere { #[typetag::serde] impl Region for Sphere { fn contains(&self, p: Meters) -> bool { - (*p - *self.center).mag() < self.rad.into_inner() + (*p - *self.center).mag() < self.rad } } @@ -140,23 +138,23 @@ impl Cube { pub fn new(lower: Meters, upper: Meters) -> Self { Self { lower, upper } } - pub fn x_range(&self) -> Range { - Flt::from_primitive(self.lower.x())..Flt::from_primitive(self.upper.x()) + pub fn x_range(&self) -> Range { + self.lower.x()..self.upper.x() } - pub fn y_range(&self) -> Range { - Flt::from_primitive(self.lower.y())..Flt::from_primitive(self.upper.y()) + pub fn y_range(&self) -> Range { + self.lower.y()..self.upper.y() } - pub fn z_range(&self) -> Range { - Flt::from_primitive(self.lower.z())..Flt::from_primitive(self.upper.z()) + pub fn z_range(&self) -> Range { + self.lower.z()..self.upper.z() } } #[typetag::serde] impl Region for Cube { fn contains(&self, p: Meters) -> bool { - self.x_range().contains(&Flt::from_primitive(p.x())) && - self.y_range().contains(&Flt::from_primitive(p.y())) && - self.z_range().contains(&Flt::from_primitive(p.z())) + self.x_range().contains(&p.x()) && + self.y_range().contains(&p.y()) && + self.z_range().contains(&p.z()) } } diff --git a/src/geom/units.rs b/src/geom/units.rs index 9966549..b5e599a 100644 --- a/src/geom/units.rs +++ b/src/geom/units.rs @@ -1,20 +1,19 @@ -use crate::flt::{Flt, Real}; -use crate::real::{Real as _, ToFloat}; +use crate::real::ToFloat; use serde::{Serialize, Deserialize}; use super::{Vec3, Vec3u}; use std::fmt::{self, Display}; use std::ops::{Add, Deref, Div, Mul, Sub}; pub trait Coord: Copy + Clone { - fn to_meters(&self, feature_size: Flt) -> Meters; - fn to_index(&self, feature_size: Flt) -> Index; - fn from_meters(other: Meters, feature_size: Flt) -> Self; - fn from_index(other: Index, feature_size: Flt) -> Self; + fn to_meters(&self, feature_size: f32) -> Meters; + fn to_index(&self, feature_size: f32) -> Index; + fn from_meters(other: Meters, feature_size: f32) -> Self; + fn from_index(other: Index, feature_size: f32) -> Self; fn from_either(i: Index, m: Meters) -> Self; } #[derive(Copy, Clone, Debug, Serialize, Deserialize)] -pub struct Meters(pub Vec3); +pub struct Meters(pub Vec3); impl Meters { pub fn new(x: R, y: R, z: R) -> Self { @@ -23,16 +22,16 @@ impl Meters { } impl Coord for Meters { - fn to_meters(&self, _feature_size: Flt) -> Meters { + fn to_meters(&self, _feature_size: f32) -> Meters { *self } - fn to_index(&self, feature_size: Flt) -> Index { - Index((self.0 / Real::from_f64(feature_size)).round().into()) + fn to_index(&self, feature_size: f32) -> Index { + Index((self.0 / feature_size).round().into()) } - fn from_meters(other: Meters, _feature_size: Flt) -> Self { + fn from_meters(other: Meters, _feature_size: f32) -> Self { other } - fn from_index(other: Index, feature_size: Flt) -> Self { + fn from_index(other: Index, feature_size: f32) -> Self { other.to_meters(feature_size) } fn from_either(_i: Index, m: Meters) -> Self { @@ -41,8 +40,8 @@ impl Coord for Meters { } impl Deref for Meters { - type Target = Vec3; - fn deref(&self) -> &Vec3 { + type Target = Vec3; + fn deref(&self) -> &Self::Target { &self.0 } } @@ -80,16 +79,16 @@ impl Index { } impl Coord for Index { - fn to_meters(&self, feature_size: Flt) -> Meters { - Meters(Vec3::from(self.0) * Real::from_f64(feature_size)) + fn to_meters(&self, feature_size: f32) -> Meters { + Meters(Vec3::from(self.0) * feature_size) } - fn to_index(&self, _feature_size: Flt) -> Index { + fn to_index(&self, _feature_size: f32) -> Index { *self } - fn from_meters(other: Meters, feature_size: Flt) -> Self { + fn from_meters(other: Meters, feature_size: f32) -> Self { other.to_index(feature_size) } - fn from_index(other: Index, _feature_size: Flt) -> Self { + fn from_index(other: Index, _feature_size: f32) -> Self { other } fn from_either(i: Index, _m: Meters) -> Self { diff --git a/src/geom/vecu.rs b/src/geom/vecu.rs index fee8746..8049c6b 100644 --- a/src/geom/vecu.rs +++ b/src/geom/vecu.rs @@ -1,5 +1,5 @@ use crate::geom::Vec3; -use crate::real::ToFloat as _; +use crate::real::Real; use serde::{Serialize, Deserialize}; use std::fmt::{self, Display}; @@ -40,8 +40,8 @@ impl From<(u32, u32, u32)> for Vec3u { } } -impl From for Vec3u { - fn from(v: Vec3) -> Self { +impl From> for Vec3u { + fn from(v: Vec3) -> Self { Self::new(v.x().to_f64() as _, v.y().to_f64() as _, v.z().to_f64() as _) } } diff --git a/src/meas.rs b/src/meas.rs index 9c7c5bc..815846b 100644 --- a/src/meas.rs +++ b/src/meas.rs @@ -1,6 +1,5 @@ -use crate::flt::{Flt, Real}; use crate::geom::{Meters, Region, Torus, Vec3, WorldRegion}; -use crate::real::Real as _; +use crate::real::{Real as _, ToFloat as _}; use crate::sim::GenericSim; use common_macros::b_tree_map; use dyn_clone::{self, DynClone}; @@ -84,19 +83,19 @@ impl Current { region: Box::new(r) } } - fn data(&self, state: &dyn GenericSim) -> (Real, Vec3) { + fn data(&self, state: &dyn GenericSim) -> (f32, Vec3) { let FieldSample(volume, current_mag, current_vec) = state.map_sum_over_enumerated(&*self.region, |coord: Meters, _cell| { let current = state.current(coord); - FieldSample(1, current.mag(), current) + FieldSample(1, current.mag().cast(), current.cast()) }); - let mean_current_mag = current_mag / (Real::from_u32(volume)); - let mean_current_vec = current_vec / (Real::from_u32(volume)); - (mean_current_mag, mean_current_vec) + let mean_current_mag = current_mag.to_f32() / (f32::from_primitive(volume)); + let mean_current_vec = current_vec.cast::() / (f32::from_primitive(volume)); + (mean_current_mag, mean_current_vec.cast()) } } #[derive(Default)] -struct FieldSample(u32, Real, Vec3); +struct FieldSample(u32, f64, Vec3); impl std::iter::Sum for FieldSample { fn sum(iter: I) -> Self @@ -172,16 +171,16 @@ impl CurrentLoop { region: r, } } - fn data(&self, state: &dyn GenericSim) -> Real { + fn data(&self, state: &dyn GenericSim) -> f32 { let FieldSample(volume, directed_current, _current_vec) = state.map_sum_over_enumerated(&self.region, |coord: Meters, _cell| { let normal = self.region.axis(); let to_coord = *coord - *self.region.center(); let tangent = normal.cross(to_coord).norm(); let current = state.current(coord); - let directed_current = current.dot(tangent); - FieldSample(1, directed_current, current) + let directed_current = current.dot(tangent.cast()); + FieldSample(1, directed_current.cast(), current.cast()) }); - let mean_directed_current = directed_current / Real::from_u32(volume); + let mean_directed_current = directed_current.cast::() / f32::from_primitive(volume); let cross_section = self.region.cross_section() / (state.feature_size() * state.feature_size()); let cross_sectional_current = mean_directed_current * cross_section; cross_sectional_current @@ -217,7 +216,7 @@ impl MagneticLoop { region: r, } } - fn data(&self, state: &dyn GenericSim) -> (Real, Real, Real) { + fn data(&self, state: &dyn GenericSim) -> (f32, f32, f32) { let FieldSamples([ FieldSample(volume, directed_m, _m_vec), FieldSample(_, directed_b, _b_vec), @@ -229,23 +228,23 @@ impl MagneticLoop { let tangent = normal.cross(to_coord).norm(); let m = cell.m(); - let directed_m = m.dot(tangent); + let directed_m = m.dot(tangent.cast()); let b = cell.b(); - let directed_b = b.dot(tangent); + let directed_b = b.dot(tangent.cast()); let h = cell.h(); - let directed_h = h.dot(tangent); + let directed_h = h.dot(tangent.cast()); FieldSamples([ - FieldSample(1, directed_m, m), - FieldSample(1, directed_b, b), - FieldSample(1, directed_h, h), + FieldSample(1, directed_m.cast(), m.cast()), + FieldSample(1, directed_b.cast(), b.cast()), + FieldSample(1, directed_h.cast(), h.cast()), ]) }); // let cross_section = self.region.cross_section() / (state.feature_size() * state.feature_size()); - let mean_directed_m = directed_m / Real::from_u32(volume); + let mean_directed_m = directed_m.cast::() / f32::from_primitive(volume); // let cross_sectional_m = mean_directed_m * cross_section; - let mean_directed_b = directed_b / Real::from_u32(volume); + let mean_directed_b = directed_b.cast::() / f32::from_primitive(volume); // let cross_sectional_b = mean_directed_b * cross_section; - let mean_directed_h = directed_h / Real::from_u32(volume); + let mean_directed_h = directed_h.cast::() / f32::from_primitive(volume); // let cross_sectional_h = mean_directed_h * cross_section; // format!( // "M({}): {:.2e}; B({}): {:.2e}; H({}): {:.2e}", @@ -292,13 +291,13 @@ impl MagneticFlux { region: Box::new(r) } } - fn data(&self, state: &dyn GenericSim) -> Vec3 { + fn data(&self, state: &dyn GenericSim) -> Vec3 { let FieldSample(volume, _directed_mag, mag_vec) = state.map_sum_over(&*self.region, |cell| { let b = cell.b(); let mag = b.mag(); - FieldSample(1, mag, b) + FieldSample(1, mag.cast(), b.cast()) }); - let mean_mag = mag_vec / Real::from_u32(volume); + let mean_mag = mag_vec.cast() / f32::from_primitive(volume); mean_mag } } @@ -331,13 +330,13 @@ impl Magnetization { region: Box::new(r) } } - fn data(&self, state: &dyn GenericSim) -> Vec3 { + fn data(&self, state: &dyn GenericSim) -> Vec3 { let FieldSample(volume, _directed_mag, mag_vec) = state.map_sum_over(&*self.region, |cell| { let m = cell.m(); let mag = m.mag(); - FieldSample(1, mag, m) + FieldSample(1, mag.cast(), m.cast()) }); - let mean_mag = mag_vec / Real::from_u32(volume); + let mean_mag = mag_vec.cast() / f32::from_primitive(volume); mean_mag } } @@ -357,7 +356,7 @@ impl AbstractMeasurement for Magnetization { } fn loc(v: Meters) -> String { - format!("{:.0} um", *v * Real::from_u64(1_000_000)) + format!("{:.0} um", *v * f32::from_primitive(1_000_000)) } /// M @@ -447,7 +446,7 @@ impl Energy { region: Box::new(region), } } - fn data(&self, state: &dyn GenericSim) -> Real { + fn data(&self, state: &dyn GenericSim) -> f32 { // Potential energy stored in a E/M field: // https://en.wikipedia.org/wiki/Magnetic_energy // https://en.wikipedia.org/wiki/Electric_potential_energy#Energy_stored_in_an_electrostatic_field_distribution @@ -456,11 +455,11 @@ impl Energy { // U(E) = 1/2 \int E . D dV #[allow(non_snake_case)] let dV = state.feature_volume(); - let e = Real::from_primitive(0.5 * dV) * state.map_sum_over(&*self.region, |cell| { + let e = f64::from_primitive(0.5 * dV) * state.map_sum_over(&*self.region, |cell| { // E . D = E . (E + P) = E.E since we don't model polarization fields - cell.h().dot(cell.b()) + cell.e().mag_sq() + cell.h().dot(cell.b()).to_f64() + cell.e().mag_sq().to_f64() }); - e + e.cast() } } @@ -494,15 +493,15 @@ impl Power { region: Box::new(region), } } - fn data(&self, state: &dyn GenericSim) -> Flt { + fn data(&self, state: &dyn GenericSim) -> f32 { // Power is P = IV = A*J*V = L^2*J.(LE) = L^3 J.E // where L is feature size. #[allow(non_snake_case)] - let dV = Real::from_primitive(state.feature_volume()); - let power = dV * state.map_sum_over(&*self.region, |cell| { - cell.current_density().dot(cell.e()) + let dV = state.feature_volume(); + let power = f64::from_primitive(dV) * state.map_sum_over(&*self.region, |cell| { + cell.current_density().dot(cell.e()).to_f64() }); - Flt::from_primitive(power) + power.cast() } } diff --git a/src/real.rs b/src/real.rs index 8c8f682..53e6a88 100644 --- a/src/real.rs +++ b/src/real.rs @@ -16,6 +16,11 @@ pub trait ToFloat { /// This exists to allow configuration over # of bits (f32 v.s. f64) as well as /// constraints. pub trait Real: ToFloat + decorum_Real + IntrinsicOrd + AddAssign + MulAssign + fmt::LowerExp + fmt::Display { + // TODO: fold with from_ + fn from_primitive(p: P) -> Self { + Self::from_f64(p.to_f64()) + } + fn from_f32(f: f32) -> Self { Self::from_f64(f as _) } @@ -23,15 +28,14 @@ pub trait Real: ToFloat + decorum_Real + IntrinsicOrd + AddAssign + MulAssign + Self::from_f32(f as _) } fn from_u64(u: u64) -> Self { - Self::from_f64(u as _) + Self::from_primitive(u) } fn from_u32(u: u32) -> Self { - Self::from_u64(u as _) + Self::from_primitive(u) } - fn from_primitive(p: P) -> Self { - // TODO: fold with from_ - Self::from_f64(p.to_f64()) + fn cast(&self) -> R { + R::from_primitive(*self) } fn tenth() -> Self { @@ -52,6 +56,23 @@ pub trait Real: ToFloat + decorum_Real + IntrinsicOrd + AddAssign + MulAssign + fn ten() -> Self { Self::from_f32(10.0) } + + fn pi() -> Self { + Self::from_f64(std::f64::consts::PI) + } + fn two_pi() -> Self { + Self::from_f64(std::f64::consts::PI * 2.0) + } + + fn c() -> Self { + Self::from_primitive(299792458) + } + fn eps0() -> Self { + Self::from_primitive(8.854187812813e-12) // F⋅m−1 + } + fn mu0() -> Self { + Self::from_primitive(1.2566370621219e-6) // H/m + } } impl ToFloat for f32 { @@ -61,12 +82,12 @@ impl ToFloat for f32 { } impl Real for f32 { + fn from_primitive(p: P) -> Self { + Self::from_f32(p.to_f32()) + } fn from_f32(f: f32) -> Self { f } - fn from_u64(u: u64) -> Self { - Self::from_f32(u as _) - } } impl ToFloat for f64 { @@ -88,12 +109,12 @@ impl ToFloat for decorum::R32 { } impl Real for decorum::R32 { + fn from_primitive(p: P) -> Self { + Self::from_f32(p.to_f32()) + } fn from_f32(f: f32) -> Self { Self::from_inner(f) } - fn from_u64(u: u64) -> Self { - Self::from_f32(u as _) - } } impl ToFloat for decorum::R64 { @@ -115,12 +136,12 @@ impl ToFloat for decorum::N32 { } impl Real for decorum::N32 { + fn from_primitive(p: P) -> Self { + Self::from_f32(p.to_f32()) + } fn from_f32(f: f32) -> Self { Self::from_inner(f) } - fn from_u64(u: u64) -> Self { - Self::from_f32(u as _) - } } impl ToFloat for decorum::N64 { @@ -135,6 +156,15 @@ impl Real for decorum::N64 { } } +impl ToFloat for i32 { + fn to_f32(&self) -> f32 { + *self as _ + } + fn to_f64(&self) -> f64 { + *self as _ + } +} + impl ToFloat for u32 { fn to_f32(&self) -> f32 { *self as _ diff --git a/src/sim.rs b/src/sim.rs index 806a992..2dec8ec 100644 --- a/src/sim.rs +++ b/src/sim.rs @@ -1,7 +1,7 @@ use crate::{flt::{Flt, Real}, consts}; 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 _}; +use crate::real::{self, Real as _, ToFloat as _}; use crate::stim::AbstractStimulus; use dyn_clone::{self, DynClone}; use log::trace; @@ -201,16 +201,16 @@ pub trait GenericSim: Send + Sync + DynClone { fn impulse_h_meters(&mut self, pos: Meters, amount: Vec3); fn impulse_b_meters(&mut self, pos: Meters, amount: Vec3); fn size(&self) -> Index; - fn feature_size(&self) -> Flt; - fn feature_volume(&self) -> Flt { + fn feature_size(&self) -> f32; + fn feature_volume(&self) -> f32 { let f = self.feature_size(); f*f*f } - fn volume(&self) -> Flt { + fn volume(&self) -> f32 { let s = self.size().to_meters(self.feature_size()); - Flt::from_primitive(s.x() * s.y() * s.z()) + s.x() * s.y() * s.z() } - fn timestep(&self) -> Flt; + fn timestep(&self) -> f32; fn step(&mut self); fn step_no(&self) -> u64; @@ -223,8 +223,8 @@ pub trait GenericSim: Send + Sync + DynClone { fn depth(&self) -> u32 { self.size().z() } - fn time(&self) -> Flt { - self.timestep() * self.step_no() as Flt + fn time(&self) -> f32 { + self.timestep() * self.step_no() as f32 } fn apply_stimulus(&mut self, stim: &dyn AbstractStimulus); @@ -389,16 +389,16 @@ impl SimState { pub fn depth(&self) -> u32 { self.size().z() } - pub fn feature_size(&self) -> Flt { - self.feature_size.into() + pub fn feature_size(&self) -> f32 { + self.feature_size.to_f32() } - pub fn timestep(&self) -> Flt { - self.timestep.into() + pub fn timestep(&self) -> f32 { + self.timestep.to_f32() } } impl SimState { - pub fn new(size: Index, feature_size: Flt) -> Self { + 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 // of a grid. The distance traveled during that time is \sqrt{3}. Hence each timestep needs @@ -406,8 +406,8 @@ impl SimState { // In 2d, this would be \sqrt{2}/2. // It's an upper limit though; it should be safe to go lower. let courant = 0.577; - let feature_size = Real::from(feature_size); - let timestep = feature_size / consts::real::C() * courant; + let feature_size = Real::from_primitive(feature_size); + let timestep = feature_size / Real::c() * Real::from_primitive(courant); Self { cells: Array3::default((size.z() as _, size.y() as _, size.x() as _)), e: Array3::default((size.z() as _, size.y() as _, size.x() as _)), @@ -456,12 +456,12 @@ impl SimState { let feature_size = self.feature_size(); let t_sec = self.time(); - let timestep = self.timestep(); + let timestep = self.timestep; trace!("apply_stimulus begin"); Zip::from(ndarray::indices_of(&self.e)).and(&mut self.e).par_apply( |(z, y, x), e| { let pos_meters = Index((x as u32, y as u32, z as u32).into()).to_meters(feature_size); - let value = stim.at(t_sec, pos_meters) * Real::from_primitive(timestep); + let value = stim.at(t_sec, pos_meters) * timestep; *e += value; }); trace!("apply_stimulus end"); @@ -507,7 +507,7 @@ impl SimState { pub fn fill_boundary_using(&mut self, thickness: C, f: F) where C: Coord, - F: Fn(Vec3) -> M, + F: Fn(Vec3) -> M, { // TODO: maybe this function belongs on the Driver? let feat = self.feature_size(); @@ -517,23 +517,23 @@ impl SimState { let region = InvertedRegion::new(Cube::new(upper_left.to_meters(feat), lower_right.to_meters(feat))); self.fill_region_using(®ion, |loc: Index| { let depth_x = if loc.x() < upper_left.x() { - (upper_left.x() - loc.x()) as Flt / upper_left.x() as Flt + (upper_left.x() - loc.x()) as f32 / upper_left.x() as f32 } else if loc.x() > lower_right.x() { - (loc.x() - lower_right.x()) as Flt / upper_left.x() as Flt + (loc.x() - lower_right.x()) as f32 / upper_left.x() as f32 } else { 0.0 }; let depth_y = if loc.y() < upper_left.y() { - (upper_left.y() - loc.y()) as Flt / upper_left.y() as Flt + (upper_left.y() - loc.y()) as f32 / upper_left.y() as f32 } else if loc.y() > lower_right.y() { - (loc.y() - lower_right.y()) as Flt / upper_left.y() as Flt + (loc.y() - lower_right.y()) as f32 / upper_left.y() as f32 } else { 0.0 }; let depth_z = if loc.z() < upper_left.z() { - (upper_left.z() - loc.z()) as Flt / upper_left.z() as Flt + (upper_left.z() - loc.z()) as f32 / upper_left.z() as f32 } else if loc.z() > lower_right.z() { - (loc.z() - lower_right.z()) as Flt / upper_left.z() as Flt + (loc.z() - lower_right.z()) as f32 / upper_left.z() as f32 } else { 0.0 }; @@ -584,11 +584,11 @@ impl GenericSim for SimState { self.cells.shape()[0] as _, )) } - fn feature_size(&self) -> Flt { - self.feature_size.into() + fn feature_size(&self) -> f32 { + self.feature_size.to_f32() } - fn timestep(&self) -> Flt { - self.timestep() + fn timestep(&self) -> f32 { + self.timestep.to_f32() } fn step(&mut self) { SimState::step(self) @@ -601,34 +601,34 @@ impl GenericSim for SimState { #[allow(unused)] impl SimState { fn get_mat(&self, c: C) -> &M { - let at = c.to_index(self.feature_size.into()); + let at = c.to_index(self.feature_size.cast()); &self.cells[at.row_major_idx()] } fn get_mat_mut(&mut self, c: C) -> &mut M { - let at = c.to_index(self.feature_size.into()); + let at = c.to_index(self.feature_size.cast()); &mut self.cells[at.row_major_idx()] } fn get_e(&self, c: C) -> &Vec3 { - let at = c.to_index(self.feature_size.into()); + let at = c.to_index(self.feature_size.cast()); &self.e[at.row_major_idx()] } fn get_e_mut(&mut self, c: C) -> &mut Vec3 { - let at = c.to_index(self.feature_size.into()); + let at = c.to_index(self.feature_size.cast()); &mut self.e[at.row_major_idx()] } fn get_h(&self, c: C) -> &Vec3 { - let at = c.to_index(self.feature_size.into()); + let at = c.to_index(self.feature_size.cast()); &self.h[at.row_major_idx()] } fn get_h_mut(&mut self, c: C) -> &mut Vec3 { - let at = c.to_index(self.feature_size.into()); + let at = c.to_index(self.feature_size.cast()); &mut self.h[at.row_major_idx()] } fn get_hcell<'a, C: Coord>(&'a mut self, c: C) -> HCell<'a, M> { - let at = c.to_index(self.feature_size.into()); + let at = c.to_index(self.feature_size.cast()); let idx = at.row_major_idx(); HCell { e: &self.e[idx], @@ -1487,14 +1487,14 @@ mod test { assert_vec3_eq(actual_df_dt, df_dt(tm)); } - fn energy(s: &dyn GenericSim) -> Flt { + fn energy(s: &dyn GenericSim) -> f32 { 0.5 * s.feature_volume() * s.map_sum_enumerated(|_pos: Index, cell| { // println!("{}: {}", _pos, cell.e()); - Flt::from_primitive(cell.e().mag_sq()) - }) + cell.e().mag_sq().to_f64() + }).to_f32() } - fn energy_now_and_then(state: &mut S, frames: u32) -> (Flt, Flt) { + fn energy_now_and_then(state: &mut S, frames: u32) -> (f32, f32) { let energy_0 = energy(state); for _f in 0..frames { // println!("e({}) = {}", _f, energy(state)); @@ -1521,7 +1521,7 @@ mod test { fn energy_conservation_over_time() { let mut state = StaticSim::new(Index((2001, 1, 1).into()), 1e-6); for t in 0..100 { - let angle = (t as Flt)/100.0*2.0*consts::PI; + let angle = (t as Flt)/100.0*Flt::two_pi(); let sig = angle.sin(); let gate = 0.5*(1.0 - angle.cos()); //println!("{}", g.exp()); @@ -1539,7 +1539,7 @@ mod test { fn sane_boundary_conditions() { let mut state = StaticSim::new(Index((21, 21, 21).into()), 1e-6); for t in 0..10 { - let angle = (t as Flt)/10.0*2.0*consts::PI; + let angle = (t as Flt)/10.0*Flt::two_pi(); let sig = angle.sin(); let gate = 0.5*(1.0 - angle.cos()); let dyn_state = &mut state as &mut dyn GenericSim; @@ -1554,11 +1554,11 @@ 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) -> (Flt, Flt) { + fn conductor_test>(mat: M) -> (f32, f32) { let mut state = StaticSim::new(Index((201, 1, 1).into()), 1e-6); state.fill_region(&WorldRegion, mat.into()); for t in 0..100 { - let angle = (t as Flt)/100.0*2.0*consts::PI; + let angle = (t as Flt)/100.0*Flt::two_pi(); let sig = angle.sin(); let gate = 0.5*(1.0 - angle.cos()); let dyn_state = &mut state as &mut dyn GenericSim; @@ -1605,8 +1605,8 @@ mod test { 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(Real::three()); - let conductivity = b * (Real::half() / timestep); + let b = boundary_ness.elem_pow(3.0); + let conductivity = b * (0.5 / timestep); // K scales up to 11; dc_shift scales from 0.05: // https://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.461.4606&rep=rep1&type=pdf // let propagation = Vec3::unit() + b * 10.0; @@ -1614,7 +1614,7 @@ mod test { // // let propagation = Vec3::unit(); // let dc_shift = Vec3::unit() * 1e-6; - Pml::new(conductivity) + Pml::new(conductivity.cast()) }); state } @@ -1624,8 +1624,8 @@ mod test { let mut state = StaticSim::new(size, 1e-6); let timestep = state.timestep(); state.fill_boundary_using(size/4, |boundary_ness| { - let b = boundary_ness.elem_pow(Real::three()); - let conductivity = Vec3::unit() * (consts::real::EPS0() * b.mag() * Real::half() / timestep); + let b = boundary_ness.elem_pow(3.0); + let conductivity = Vec3::uniform(f32::eps0() * b.mag() * 0.5 / timestep); Static { conductivity, ..Default::default() } @@ -1635,14 +1635,15 @@ mod test { /// Stimulus which applies a (possibly different) e*sin(wt)*gate(t) at every point struct PmlStim { + /// Maps location -> (field direction, frequency) f: F, - t_end: Flt, - feat_size: Flt, - sqrt_vol: Flt, + t_end: f32, + feat_size: f32, + sqrt_vol: f32, } - impl (Vec3, Flt) + Sync> AbstractStimulus for PmlStim { - fn at(&self, t_sec: Flt, pos: Meters) -> Vec3 { - let angle = (t_sec as Flt)/self.t_end*2.0*consts::PI; + impl (Vec3, f32) + Sync> AbstractStimulus for PmlStim { + fn at(&self, t_sec: f32, pos: Meters) -> Vec3 { + let angle = t_sec/self.t_end * f32::two_pi(); let gate = 0.5*(1.0 - angle.cos()); let (e, hz) = (self.f)(pos.to_index(self.feat_size)); let sig_angle = angle*hz; @@ -1652,8 +1653,8 @@ mod test { } /// Returns the energy ratio (Eend / Estart) - fn pml_test_full_interior(state: &mut S, f: F) -> Flt - where F: Fn(Index) -> (Vec3, Flt) + Sync // returns (E vector, cycles (like Hz)) + fn pml_test_full_interior(state: &mut S, f: F) -> f32 + where F: Fn(Index) -> (Vec3, f32) + Sync // returns (E vector, cycles (like Hz)) { let stim = PmlStim { f, @@ -1670,13 +1671,13 @@ mod test { // sanity check the starting energy assert_gt!(energy_0, 1e-11); assert_lt!(energy_0, 1.0); - (Real::from(energy_1) / Real::from(energy_0)).into_inner() + energy_1 / energy_0 } - fn pml_test_over_region(state: &mut S, reg: R, f: F) -> Flt + fn pml_test_over_region(state: &mut S, reg: R, f: F) -> f32 where R: Region, - F: Fn(Index) -> (Vec3, Flt) + Sync, + F: Fn(Index) -> (Vec3, f32) + Sync, S: GenericSim { let feat = state.feature_size(); @@ -1688,13 +1689,13 @@ mod test { } }) } - fn pml_test_uniform_over_region(state: &mut S, reg: R, e: Vec3, hz: Flt) -> Flt { + fn pml_test_uniform_over_region(state: &mut S, reg: R, e: Vec3, hz: f32) -> f32 { pml_test_over_region(state, reg, |_idx| { (e, hz) }) } - fn pml_test_at(state: &mut S, e: Vec3, center: Index) -> Flt { + fn pml_test_at(state: &mut S, e: Vec3, center: Index) -> f32 { pml_test_full_interior(state, |idx| { if idx == center { (e, 1.0) @@ -1770,17 +1771,17 @@ mod test { 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(Real::three()); - let coord_stretch = b * (Real::half() / timestep); + let b = boundary_ness.elem_pow(3.0); + let coord_stretch = b * (0.5 / timestep); // let coord_stretch = b * 0.5; // Force the stretch to be only along one axis - let coord_stretch = match coord_stretch.to_tuple() { - (x, y, z) if x >= y && x >= z => (x, Real::zero(), Real::zero()), - (x, y, z) if y >= x && y >= z => (Real::zero(), y, Real::zero()), - (x, y, z) if z >= x && z >= y => (Real::zero(), Real::zero(), z), + let coord_stretch: Vec3<_> = match coord_stretch.to_tuple() { + (x, y, z) if x >= y && x >= z => (x, 0.0, 0.0), + (x, y, z) if y >= x && y >= z => (0.0, y, 0.0), + (x, y, z) if z >= x && z >= y => (0.0, 0.0, z), _ => unreachable!(), }.into(); - Pml::new(coord_stretch) + Pml::new(coord_stretch.cast()) }); state } @@ -1792,11 +1793,11 @@ mod test { 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(Real::three()); - let cs = b * (Real::half() / timestep); + let b = boundary_ness.elem_pow(3.0); + let cs = b * (0.5 / timestep); // let cs = b * 0.5; // permute the stretching so that it shouldn't interfere with the wave - let coord_stretch = shuffle(cs); + let coord_stretch = shuffle(cs.cast()); Pml::new(coord_stretch) }); let center = Index(state.size().0 / 2); diff --git a/src/stim.rs b/src/stim.rs index 282cbc2..73b4cd9 100644 --- a/src/stim.rs +++ b/src/stim.rs @@ -1,27 +1,26 @@ -use crate::consts; use crate::flt::{self, Flt}; use crate::real::*; use crate::geom::{Meters, Region, Vec3}; pub trait AbstractStimulus: Sync { /// Return the E field which should be added PER-SECOND to the provided position/time. - fn at(&self, t_sec: Flt, pos: Meters) -> Vec3; + fn at(&self, t_sec: f32, pos: Meters) -> Vec3; } // impl AbstractStimulus for &T { -// fn at(&self, t_sec: Flt, pos: Meters) -> Vec3 { +// fn at(&self, t_sec: f32, pos: Meters) -> Vec3 { // (*self).at(t_sec, pos) // } // } impl AbstractStimulus for Vec { - fn at(&self, t_sec: Flt, pos: Meters) -> Vec3 { + fn at(&self, t_sec: f32, pos: Meters) -> Vec3 { self.iter().map(|s| s.at(t_sec, pos)).sum() } } impl AbstractStimulus for Box { - fn at(&self, t_sec: Flt, pos: Meters) -> Vec3 { + fn at(&self, t_sec: f32, pos: Meters) -> Vec3 { (**self).at(t_sec, pos) } } @@ -42,7 +41,7 @@ impl Stimulus { } impl AbstractStimulus for Stimulus { - fn at(&self, t_sec: Flt, pos: Meters) -> Vec3 { + fn at(&self, t_sec: f32, pos: Meters) -> Vec3 { if self.region.contains(pos) { self.stim.at(t_sec) } else { @@ -68,12 +67,12 @@ impl CurlStimulus { } impl AbstractStimulus for CurlStimulus { - fn at(&self, t_sec: Flt, pos: Meters) -> Vec3 { + fn at(&self, t_sec: f32, pos: Meters) -> Vec3 { if self.region.contains(pos) { let amount = self.stim.at(t_sec); let from_center_to_point = *pos - *self.center; - let rotational = from_center_to_point.cross(*self.axis); - let impulse = rotational.with_mag(flt::Real::from_primitive(amount)); + let rotational: Vec3 = from_center_to_point.cross(*self.axis).cast(); + let impulse = rotational.with_mag(amount.cast()); impulse } else { Vec3::zero() @@ -84,22 +83,22 @@ impl AbstractStimulus for CurlStimulus pub trait TimeVarying1: Sized { /// Retrieve the E impulse to apply PER-SECOND at the provided time (in seconds). - fn at(&self, t_sec: Flt) -> Flt; - fn shifted(self, new_start: Flt) -> Shifted { + fn at(&self, t_sec: f32) -> Flt; + fn shifted(self, new_start: f32) -> Shifted { Shifted::new(self, new_start) } - fn gated(self, from: Flt, to: Flt) -> Gated { + fn gated(self, from: f32, to: f32) -> Gated { Gated::new(self, from, to) } } pub trait TimeVarying3: Sized { /// Retrieve the E impulse to apply PER-SECOND at the provided time (in seconds). - fn at(&self, t_sec: Flt) -> Vec3; - fn shifted(self, new_start: Flt) -> Shifted { + fn at(&self, t_sec: f32) -> Vec3; + fn shifted(self, new_start: f32) -> Shifted { Shifted::new(self, new_start) } - fn gated(self, from: Flt, to: Flt) -> Gated { + fn gated(self, from: f32, to: f32) -> Gated { Gated::new(self, from, to) } } @@ -107,26 +106,26 @@ pub trait TimeVarying3: Sized { #[derive(Clone)] pub struct Sinusoid { amp: A, - omega: Flt, + omega: f32, } pub type Sinusoid1 = Sinusoid; pub type Sinusoid3 = Sinusoid; impl Sinusoid { - pub fn new(amp: A, freq: Flt) -> Self { + pub fn new(amp: A, freq: f32) -> Self { Self { amp, - omega: freq * consts::TWO_PI, + omega: freq * f32::two_pi(), } } - pub fn from_wavelength(amp: A, lambda: Flt) -> Self { + pub fn from_wavelength(amp: A, lambda: f32) -> Self { Self::new(amp, 1.0/lambda) } - pub fn freq(&self) -> Flt { - self.omega / consts::TWO_PI + pub fn freq(&self) -> f32 { + self.omega / f32::two_pi() } - pub fn wavelength(&self) -> Flt { + pub fn wavelength(&self) -> f32 { 1.0 / self.freq() } pub fn one_cycle(self) -> Gated { @@ -140,13 +139,13 @@ impl Sinusoid { } impl TimeVarying1 for Sinusoid1 { - fn at(&self, t_sec: Flt) -> Flt { - self.amp * (t_sec * self.omega).sin() + fn at(&self, t_sec: f32) -> Flt { + self.amp * (t_sec * self.omega).cast::().sin() } } impl TimeVarying3 for Sinusoid3 { - fn at(&self, t_sec: Flt) -> Vec3 { + fn at(&self, t_sec: f32) -> Vec3 { self.amp * flt::Real::from_primitive(t_sec * self.omega).sin() } } @@ -154,18 +153,18 @@ impl TimeVarying3 for Sinusoid3 { #[derive(Clone)] pub struct Gated { inner: T, - start: Flt, - end: Flt, + start: f32, + end: f32, } impl Gated { - pub fn new(inner: T, start: Flt, end: Flt) -> Self { + pub fn new(inner: T, start: f32, end: f32) -> Self { Self { inner, start, end } } } impl TimeVarying1 for Gated { - fn at(&self, t_sec: Flt) -> Flt { + fn at(&self, t_sec: f32) -> Flt { if (self.start..self.end).contains(&t_sec) { self.inner.at(t_sec) } else { @@ -175,7 +174,7 @@ impl TimeVarying1 for Gated { } impl TimeVarying3 for Gated { - fn at(&self, t_sec: Flt) -> Vec3 { + fn at(&self, t_sec: f32) -> Vec3 { if (self.start..self.end).contains(&t_sec) { self.inner.at(t_sec) } else { @@ -187,23 +186,23 @@ impl TimeVarying3 for Gated { #[derive(Clone)] pub struct Shifted { inner: T, - start_at: Flt, + start_at: f32, } impl Shifted { - pub fn new(inner: T, start_at: Flt) -> Self { + pub fn new(inner: T, start_at: f32) -> Self { Self { inner, start_at } } } impl TimeVarying1 for Shifted { - fn at(&self, t_sec: Flt) -> Flt { + fn at(&self, t_sec: f32) -> Flt { self.inner.at(t_sec - self.start_at) } } impl TimeVarying3 for Shifted { - fn at(&self, t_sec: Flt) -> Vec3 { + fn at(&self, t_sec: f32) -> Vec3 { self.inner.at(t_sec - self.start_at) } }