From dfdbb4318049e68f282a007f6ccfc8fb44f3c601 Mon Sep 17 00:00:00 2001 From: Colin Date: Tue, 8 Jun 2021 18:25:10 -0700 Subject: [PATCH] Parameterize most of SimState and Material over Real. Some things still work only for coremem::Real. Need to troubleshoot those. --- benches/driver.rs | 9 +- examples/minimal_torus.rs | 4 +- examples/toroid25d.rs | 4 +- examples/wrapped_torus.rs | 4 +- src/bin/diagnostics.rs | 2 +- src/driver.rs | 6 +- src/geom/vec.rs | 2 +- src/mat.rs | 430 +++++++++++++++++++++----------------- src/post.rs | 2 +- src/real.rs | 2 +- src/sim.rs | 64 +++--- 11 files changed, 285 insertions(+), 244 deletions(-) diff --git a/benches/driver.rs b/benches/driver.rs index ae32b8a..9ede69b 100644 --- a/benches/driver.rs +++ b/benches/driver.rs @@ -1,10 +1,11 @@ +use coremem::Real; use coremem::{Driver, geom::Index, mat::GenericMaterial, mat::GenericMaterialNoPml, mat::GenericMaterialOneField}; use criterion::{BenchmarkId, criterion_group, criterion_main, Criterion}; pub fn bench_step(c: &mut Criterion) { for size in &[10, 20, 40, 80, 160] { c.bench_with_input(BenchmarkId::new("Driver::step", size), size, |b, &size| { - let mut driver = Driver::::new(Index::new(size, size, size), 1e-5); + let mut driver = Driver::>::new(Index::new(size, size, size), 1e-5); b.iter(|| driver.step()) }); } @@ -13,7 +14,7 @@ pub fn bench_step(c: &mut Criterion) { pub fn bench_step_no_pml(c: &mut Criterion) { for size in &[10, 20, 40, 80, 160] { c.bench_with_input(BenchmarkId::new("Driver::step_no_pml", size), size, |b, &size| { - let mut driver = Driver::::new(Index::new(size, size, size), 1e-5); + let mut driver = Driver::>::new(Index::new(size, size, size), 1e-5); b.iter(|| driver.step()) }); } @@ -22,7 +23,7 @@ pub fn bench_step_no_pml(c: &mut Criterion) { pub fn bench_step_one_vec(c: &mut Criterion) { for size in &[10, 20, 40, 80, 160] { c.bench_with_input(BenchmarkId::new("Driver::step_one_vec", size), size, |b, &size| { - let mut driver = Driver::::new(Index::new(size, size, size), 1e-5); + let mut driver = Driver::>::new(Index::new(size, size, size), 1e-5); b.iter(|| driver.step()) }); } @@ -32,7 +33,7 @@ pub fn bench_step_with_pml(c: &mut Criterion) { let size = 40; for thickness in &[0, 1, 2, 4, 8, 16] { c.bench_with_input(BenchmarkId::new("Driver::step_with_pml", thickness), thickness, |b, &thickness| { - let mut driver = Driver::::new(Index::new(size, size, size), 1e-5); + let mut driver = Driver::>::new(Index::new(size, size, size), 1e-5); driver.add_pml_boundary(Index::new(thickness, thickness, thickness)); b.iter(|| driver.step()) }); diff --git a/examples/minimal_torus.rs b/examples/minimal_torus.rs index bc9b833..4a52958 100644 --- a/examples/minimal_torus.rs +++ b/examples/minimal_torus.rs @@ -1,4 +1,4 @@ -use coremem::{Driver, Flt, mat, meas}; +use coremem::{Driver, Flt, Real, mat, meas}; use coremem::geom::{Cube, Index, InvertedRegion, Meters, Torus, Union}; use coremem::stim::{CurlStimulus, Sinusoid1, TimeVarying1 as _}; @@ -26,7 +26,7 @@ fn main() { let width_px = from_m(width); let depth_px = from_m(depth); let size_px = Index((width_px, width_px, depth_px).into()); - let mut driver: Driver = Driver::new(size_px, feat_size); + let mut driver: Driver> = Driver::new(size_px, feat_size); driver.set_steps_per_frame(400); driver.set_steps_per_stim(1); let base = "minimal_torus-6"; diff --git a/examples/toroid25d.rs b/examples/toroid25d.rs index ce2951f..e8d8ae2 100644 --- a/examples/toroid25d.rs +++ b/examples/toroid25d.rs @@ -1,4 +1,4 @@ -use coremem::{Driver, Flt, mat, meas}; +use coremem::{Driver, Flt, Real, mat, meas}; use coremem::geom::{CylinderZ, Index, Meters, Vec2, Vec3}; use coremem::stim::{Stimulus, Sinusoid3}; use log::trace; @@ -30,7 +30,7 @@ fn main() { let width_px = from_m(width); let depth_px = from_m(depth); let size_px = Index((width_px, width_px, depth_px).into()); - let mut driver: Driver = Driver::new(size_px, feat_size); + let mut driver: Driver> = Driver::new(size_px, feat_size); //driver.set_steps_per_frame(8); //driver.set_steps_per_frame(20); //driver.set_steps_per_frame(40); diff --git a/examples/wrapped_torus.rs b/examples/wrapped_torus.rs index d083531..9c81e74 100644 --- a/examples/wrapped_torus.rs +++ b/examples/wrapped_torus.rs @@ -1,4 +1,4 @@ -use coremem::{Driver, Flt, mat, meas}; +use coremem::{Driver, Flt, Real, mat, meas}; use coremem::geom::{Index, Meters, Torus}; use coremem::stim::{CurlStimulus, Sinusoid1, TimeVarying1 as _}; @@ -32,7 +32,7 @@ fn main() { let height_px = from_m(height); let depth_px = from_m(depth); let size_px = Index((width_px, height_px, depth_px).into()); - let mut driver: Driver = Driver::new(size_px, feat_size); + let mut driver: Driver> = Driver::new(size_px, feat_size); driver.set_steps_per_frame(500); // driver.set_steps_per_stim(10); let base = "wrapped_torus-21-high-input-resistance"; diff --git a/src/bin/diagnostics.rs b/src/bin/diagnostics.rs index 66b640e..5c31e33 100644 --- a/src/bin/diagnostics.rs +++ b/src/bin/diagnostics.rs @@ -1,7 +1,7 @@ use coremem::mat::*; fn main() { - let m3r1 = Ferroxcube3R1::curve(); + let m3r1 = Ferroxcube3R1::::curve(); let extremes = m3r1.extremes(); println!("3R1 Extremes: H={}, M={}", extremes.x(), extremes.y()); } diff --git a/src/driver.rs b/src/driver.rs index fbdc653..5964367 100644 --- a/src/driver.rs +++ b/src/driver.rs @@ -15,7 +15,7 @@ use std::sync::mpsc::{sync_channel, SyncSender, Receiver}; use std::time::{Duration, Instant}; use threadpool::ThreadPool; -pub struct Driver { +pub struct Driver, R=flt::Real> { pub state: SimState, renderer: Arc>>, // TODO: use Rayon's thread pool? @@ -192,7 +192,7 @@ impl + Clone + Default + Send + } } -impl + From> Driver { +impl + From>> Driver { pub fn add_pml_boundary(&mut self, thickness: C) { let timestep = self.state.timestep(); self.state.fill_boundary_using(thickness, |boundary_ness| { @@ -203,7 +203,7 @@ impl + From> Driver { } } -impl + From> Driver { +impl + From>> Driver { pub fn add_classical_boundary(&mut self, thickness: C) { let timestep = self.state.timestep(); self.state.fill_boundary_using(thickness, |boundary_ness| { diff --git a/src/geom/vec.rs b/src/geom/vec.rs index 84c3213..ae037ee 100644 --- a/src/geom/vec.rs +++ b/src/geom/vec.rs @@ -7,7 +7,7 @@ use std::fmt; use std::iter::Sum; use std::ops::{Add, AddAssign, Div, DivAssign, Mul, MulAssign, Neg, Sub}; -#[derive(Copy, Clone, Debug, Default, Serialize, Deserialize)] +#[derive(Copy, Clone, Debug, Default, Eq, PartialEq, Serialize, Deserialize)] pub struct Vec2 { pub x: R, pub y: R, diff --git a/src/mat.rs b/src/mat.rs index 1c7c5a3..6c4bcb3 100644 --- a/src/mat.rs +++ b/src/mat.rs @@ -8,10 +8,13 @@ use lazy_static::lazy_static; use log::trace; use enum_dispatch::enum_dispatch; use serde::{Serialize, Deserialize}; +use std::any::{Any, TypeId}; use std::cmp::Ordering; +use std::collections::HashMap; +use std::sync::Mutex; #[enum_dispatch] -pub trait Material { +pub trait Material { fn step_parameters_mut<'a>(&'a mut self) -> StepParametersMut<'a, R> { // by default, behave as a vacuum StepParametersMut::default() @@ -42,14 +45,14 @@ impl> MaterialExt for M { /// Capable of capturing all field-related information about a material at any /// snapshot moment-in-time. Useful for serializing state. #[derive(Clone, Default, Serialize, Deserialize)] -pub struct Static { - pub conductivity: Vec3, +pub struct Static { + pub conductivity: Vec3, // pub pml: Option<(PmlState, PmlParameters)>, - pub m: Vec3, + pub m: Vec3, } -impl Static { - pub fn from_material>(m: &M) -> Self { +impl Static { + pub fn from_material>(m: &M) -> Self { let p = m.step_parameters(); Self { conductivity: p.conductivity(), @@ -63,20 +66,20 @@ impl Static { // } } -impl Material for Static { - fn step_parameters_mut<'a>(&'a mut self) -> StepParametersMut<'a, flt::Real> { +impl Material for Static { + fn step_parameters_mut<'a>(&'a mut self) -> StepParametersMut<'a, R> { StepParametersMut::new( self.conductivity, None, // self.pml.as_mut().map(|(s, p)| (s, *p)), ) } - fn m(&self) -> Vec3 { + fn m(&self) -> Vec3 { self.m } } -impl From for Static -where T: Into +impl From for Static +where T: Into> { fn from(mat: T) -> Self { let generic = mat.into(); @@ -86,45 +89,45 @@ where T: Into /// Material which has a conductivity parameter, but cannot be magnetized #[derive(Clone, Default, Serialize, Deserialize)] -pub struct Conductor { - conductivity: Vec3, +pub struct Conductor { + conductivity: Vec3, } -impl Conductor { - pub fn new(conductivity: R) -> Self { +impl Conductor { + pub fn new(conductivity: R2) -> Self { Self { conductivity: Vec3::uniform(conductivity).cast() } } - pub fn new_anisotropic(conductivity: Vec3) -> Self { + pub fn new_anisotropic(conductivity: Vec3) -> Self { Self { conductivity: conductivity.cast(), } } } -impl Material for Conductor { - fn step_parameters_mut<'a>(&'a mut self) -> StepParametersMut<'a, flt::Real> { +impl Material for Conductor { + fn step_parameters_mut<'a>(&'a mut self) -> StepParametersMut<'a, R> { StepParametersMut::default().with_conductivity(self.conductivity) } } /// Material which can be magnetized, but has no hysteresis and no coercivity. #[derive(Clone, Default, Serialize, Deserialize)] -pub struct LinearMagnet { +pub struct LinearMagnet { /// \mu_r - relative_permeability: Vec3, - m: Vec3, + relative_permeability: Vec3, + m: Vec3, } -impl LinearMagnet { - pub fn new(relative_permeability: R) -> Self { +impl LinearMagnet { + pub fn new(relative_permeability: R2) -> Self { Self { relative_permeability: Vec3::uniform(relative_permeability).cast(), m: Vec3::zero(), } } - pub fn new_anisotropic(relative_permeability: Vec3) -> Self { + pub fn new_anisotropic(relative_permeability: Vec3) -> Self { Self { relative_permeability: relative_permeability.cast(), m: Vec3::zero() @@ -132,11 +135,11 @@ impl LinearMagnet { } } -impl Material for LinearMagnet { - fn m(&self) -> Vec3 { +impl Material for LinearMagnet { + fn m(&self) -> Vec3 { self.m } - fn step_b(&mut self, _context: &CellState, delta_b: Vec3) { + fn step_b(&mut self, _context: &CellState, delta_b: Vec3) { //```tex // $B = \mu_0 (H + M) = \mu_0 \mu_r H$ // $\mu_r H = H + M$ @@ -145,80 +148,107 @@ impl Material for LinearMagnet { // $B = \mu_0 \mu_r/(\mu_r - 1) M$ //``` let mu_r = self.relative_permeability; - let delta_m = (delta_b*flt::Real::mu0_inv()).elem_mul(mu_r - Vec3::unit()).elem_div(mu_r); + let delta_m = (delta_b*R::mu0_inv()).elem_mul(mu_r - Vec3::unit()).elem_div(mu_r); self.m += delta_m; } } #[derive(Clone, Default, Serialize, Deserialize)] -pub struct Pml(PmlState, PmlParameters); +pub struct Pml(PmlState, PmlParameters); -impl Pml { - pub fn new(pseudo_conductivity: Vec3) -> Self { +impl Pml { + pub fn new(pseudo_conductivity: Vec3) -> Self { Self(PmlState::new(), PmlParameters::new(pseudo_conductivity)) } } -impl Material for Pml { - fn step_parameters_mut<'a>(&'a mut self) -> StepParametersMut<'a, flt::Real> { +impl Material for Pml { + fn step_parameters_mut<'a>(&'a mut self) -> StepParametersMut<'a, R> { StepParametersMut::default().with_pml(&mut self.0, self.1) } } -pub trait PiecewiseLinearFerromagnet { - fn curve() -> &'static MHCurve; - fn conductivity() -> Flt; - fn m(&self) -> Vec3; - fn m_mut(&mut self) -> &mut Vec3; -} +// pub trait PiecewiseLinearFerromagnet { +// fn curve() -> &'static MHCurve; +// fn conductivity() -> R; +// fn m(&self) -> Vec3; +// fn m_mut(&mut self) -> &mut Vec3; +// } +// +// impl> Material for T { +// fn step_parameters_mut<'a>(&'a mut self) -> StepParametersMut<'a, R> { +// let c = T::conductivity(); +// StepParametersMut::default().with_conductivity(Vec3::uniform(c)) +// } +// fn m(&self) -> Vec3 { +// self.m() +// } +// fn step_b(&mut self, context: &CellState, delta_b: Vec3) { +// trace!("step_b enter"); +// let mh_curve = T::curve(); +// let (h, m) = (context.h(), self.m()); +// let target_hm = h + m + delta_b * R::mu0_inv(); +// +// // TODO: this is probably not the best way to generalize a BH curve into 3d. +// let (_hx, mx) = mh_curve.move_to( +// h.x(), +// m.x(), +// target_hm.x(), +// ); +// let (_hy, my) = mh_curve.move_to( +// h.y(), +// m.y(), +// target_hm.y(), +// ); +// let (_hz, mz) = mh_curve.move_to( +// h.z(), +// m.z(), +// target_hm.z(), +// ); +// +// *self.m_mut() = Vec3::new(mx, my, mz); +// // let ret = Vec3::new(hx, hy, hz); +// trace!("step_b end"); +// } +// } -impl Material for T { - fn step_parameters_mut<'a>(&'a mut self) -> StepParametersMut<'a, flt::Real> { - let c = T::conductivity(); - StepParametersMut::default().with_conductivity(Vec3::uniform(c)) - } - fn m(&self) -> Vec3 { - self.m() - } - fn step_b(&mut self, context: &CellState, delta_b: Vec3) { - trace!("step_b enter"); - let mh_curve = T::curve(); - let (h, m) = (context.h(), self.m()); - let target_hm = h + m + delta_b * flt::Real::mu0_inv(); +fn step_linear_ferro(m_mut: &mut Vec3, mh_curve: &MHCurve, context: &CellState, delta_b: Vec3) { + trace!("step_b enter"); + let (h, m) = (context.h(), *m_mut); + let target_hm = h + m + delta_b * R::mu0_inv(); - // TODO: this is probably not the best way to generalize a BH curve into 3d. - let (_hx, mx) = mh_curve.move_to( - Flt::from_primitive(h.x()), - Flt::from_primitive(m.x()), - Flt::from_primitive(target_hm.x()) - ); - let (_hy, my) = mh_curve.move_to( - Flt::from_primitive(h.y()), - Flt::from_primitive(m.y()), - Flt::from_primitive(target_hm.y()) - ); - let (_hz, mz) = mh_curve.move_to( - Flt::from_primitive(h.z()), - Flt::from_primitive(m.z()), - Flt::from_primitive(target_hm.z()) - ); + // TODO: this is probably not the best way to generalize a BH curve into 3d. + let (_hx, mx) = mh_curve.move_to( + h.x(), + m.x(), + target_hm.x(), + ); + let (_hy, my) = mh_curve.move_to( + h.y(), + m.y(), + target_hm.y(), + ); + let (_hz, mz) = mh_curve.move_to( + h.z(), + m.z(), + target_hm.z(), + ); - *self.m_mut() = Vec3::new(mx, my, mz).cast(); - // let ret = Vec3::new(hx, hy, hz); - trace!("step_b end"); - } + *m_mut = Vec3::new(mx, my, mz); + // let ret = Vec3::new(hx, hy, hz); + trace!("step_b end"); } /// M as a function of H #[derive(Clone)] -pub struct MHCurve { - geom: Polygon2d, +pub struct MHCurve { + geom: Polygon2d, } -impl MHCurve { +impl MHCurve { /// Construct a M(H) curve from a sweep from M = 0 to Ms and back down to M = 0. /// The curve below M = 0 is derived by symmetry. - fn new(points: &[Vec2]) -> Self { + fn new(points: &[Vec2]) -> Self { let full_pts: Vec<_> = points.iter().cloned() .chain(points.iter().cloned().map(|p| -p)) @@ -230,15 +260,15 @@ impl MHCurve { } } - fn from_bh(points: &[(R, R)]) -> Self { + fn from_bh(points: &[(R2, R2)]) -> Self { let mh_points: Vec<_> = points.iter().cloned().map(|(h, b)| { - Vec2::new(h, b / R::mu0() - h) + Vec2::new(h, b / R2::mu0() - h) }).collect(); Self::new(&*mh_points) } - fn from_mh(points: &[(R, R)]) -> Self { + fn from_mh(points: &[(R2, R2)]) -> Self { let mh_points: Vec<_> = points.iter().cloned().map(|(h, m)| { Vec2::new(h, m) }).collect(); @@ -247,22 +277,22 @@ impl MHCurve { } /// Return (Hmax, Mmax) - pub fn extremes(&self) -> Vec2 { + pub fn extremes(&self) -> Vec2 { Vec2::new(self.geom.max_x(), self.geom.max_y()) } /// Moves (h, m) towards some location in the MH curve where H + M = target_hm. /// Returns `Ok((h, m))` if complete; `Err((h, m))` if there's more work to be done (call it /// again). - fn step_toward(&self, h: flt::Real, m: flt::Real, target_hm: flt::Real) -> Result<(flt::Real, flt::Real), (flt::Real, flt::Real)> { + fn step_toward(&self, h: R, m: R, target_hm: R) -> Result, Vec2> { let is_ascending = match target_hm.partial_cmp(&(h + m)).unwrap_or_else(|| panic!("{} {}", h, m)) { Ordering::Greater => true, Ordering::Less => false, - _ => return Ok((h, m)) + _ => return Ok(Vec2::new(h, m)) }; if (is_ascending && m == self.geom.max_y()) || (!is_ascending && m == self.geom.min_y()) { // Fully saturated. m is fixed, while h moves freely - return Ok((target_hm - m, m)); + return Ok(Vec2::new(target_hm - m, m)); } // Locate the segment which would contain the current point let mut segments = self.geom.segments(); @@ -271,7 +301,7 @@ impl MHCurve { panic!("failed to find segment for h:{}, m:{}, {:?}", h, m, self.geom.segments().collect::>()); }); if line.contains_y(m) && line.is_ascending() == is_ascending { - if line.contains_x(h) && line.distance_sq(Vec2::new(h, m)) < 1.0e-6 { + if line.contains_x(h) && line.distance_sq(Vec2::new(h, m)) < R::from_primitive(1.0e-6) { // (h, m) resides on this line break line; } else { @@ -298,22 +328,21 @@ impl MHCurve { if sum_h.contains_x(new_h) { // the segment contains a point with the target H+M - Ok(active_segment.at_x(new_h).into()) + Ok(active_segment.at_x(new_h)) } else { // the segment doesn't contain the desired point: clamp and try the next segment - Err(active_segment.clamp_by_x(new_h).into()) + Err(active_segment.clamp_by_x(new_h)) } } - fn move_to(&self, h: Flt, m: Flt, target_hm: Flt) -> (Flt, Flt) { + fn move_to(&self, mut h: R, mut m: R, target_hm: R) -> (R, R) { let mut i = 0; - let (mut h, mut m, target_hm) = (flt::Real::from_inner(h), flt::Real::from_inner(m), flt::Real::from_inner(target_hm)); loop { i += 1; match self.step_toward(h, m, target_hm) { - Ok((x, y)) => break (x.into(), y.into()), - Err((x, y)) => { - h = x; - m = y; + Ok(v) => break (v.x(), v.y()), + Err(v) => { + h = v.x(); + m = v.y(); }, } if i % 2048 == 0 { @@ -324,14 +353,18 @@ impl MHCurve { } #[derive(Default, Copy, Clone, Serialize, Deserialize)] -pub struct Ferroxcube3R1 { - m: Vec3, +pub struct Ferroxcube3R1 { + m: Vec3, } -impl PiecewiseLinearFerromagnet for Ferroxcube3R1 { - fn curve() -> &'static MHCurve { +impl Ferroxcube3R1 { + pub fn curve() -> &'static MHCurve { lazy_static! { - static ref FERROXCUBE_3R1: MHCurve = MHCurve::from_bh(&[ + static ref curves: Mutex>> = Mutex::new(HashMap::new()); + } + let mut lock = curves.lock().unwrap(); + let curve = lock.entry(TypeId::of::()).or_insert_with(|| { + Box::new(MHCurve::::from_bh(&[ ( 35.0, 0.0), ( 50.0, 0.250), ( 100.0, 0.325), @@ -342,99 +375,109 @@ impl PiecewiseLinearFerromagnet for Ferroxcube3R1 { ( 100.0, 0.345), ( 50.0, 0.340), ( 0.0, 0.325), - ]); - } - &*FERROXCUBE_3R1 + ])) + }).downcast_ref::>().unwrap(); + unsafe { std::mem::transmute::<&MHCurve, &'static MHCurve>(curve) } } - fn conductivity() -> Flt { - 1e-3 +} + +impl Material for Ferroxcube3R1 { + fn step_b(&mut self, context: &CellState, delta_b: Vec3) { + step_linear_ferro(&mut self.m, Self::curve(), context, delta_b) } - fn m(&self) -> Vec3 { + fn m(&self) -> Vec3 { self.m } - fn m_mut(&mut self) -> &mut Vec3 { - &mut self.m + fn step_parameters_mut<'a>(&'a mut self) -> StepParametersMut<'a, R> { + StepParametersMut::default().with_conductivity(Vec3::uniform(1e-3)) } } /// Simple, square-loop ferrite #[derive(Default, Copy, Clone, Serialize, Deserialize)] -pub struct MinimalSquare { - m: Vec3, +pub struct MinimalSquare { + m: Vec3, } -impl PiecewiseLinearFerromagnet for MinimalSquare { - fn curve() -> &'static MHCurve { +impl MinimalSquare { + pub fn curve() -> &'static MHCurve { lazy_static! { - static ref CURVE: MHCurve = MHCurve::from_mh(&[ + static ref curves: Mutex>> = Mutex::new(HashMap::new()); + } + let mut lock = curves.lock().unwrap(); + let curve = lock.entry(TypeId::of::()).or_insert_with(|| { + Box::new(MHCurve::::from_bh(&[ ( 1.0, 0.0), ( 2.0, 1000000.0), // Falling ( 0.0, 900000.0), - ]); - } - &*CURVE + ])) + }).downcast_ref::>().unwrap(); + unsafe { std::mem::transmute::<&MHCurve, &'static MHCurve>(curve) } } - fn conductivity() -> Flt { - 1e-3 +} + +impl Material for MinimalSquare { + fn step_b(&mut self, context: &CellState, delta_b: Vec3) { + step_linear_ferro(&mut self.m, Self::curve(), context, delta_b) } - fn m(&self) -> Vec3 { + fn m(&self) -> Vec3 { self.m } - fn m_mut(&mut self) -> &mut Vec3 { - &mut self.m + fn step_parameters_mut<'a>(&'a mut self) -> StepParametersMut<'a, R> { + StepParametersMut::default().with_conductivity(Vec3::uniform(1e-3)) } } // #[enum_dispatch(Material)] #[derive(Clone, Serialize, Deserialize)] -pub enum GenericMaterial { - Conductor(Conductor), - LinearMagnet(LinearMagnet), - Pml(Pml), - Ferroxcube3R1(Ferroxcube3R1), - MinimalSquare(MinimalSquare), +pub enum GenericMaterial { + Conductor(Conductor), + LinearMagnet(LinearMagnet), + Pml(Pml), + Ferroxcube3R1(Ferroxcube3R1), + MinimalSquare(MinimalSquare), } -impl Default for GenericMaterial { +impl Default for GenericMaterial { fn default() -> Self { Conductor::default().into() } } -impl From for GenericMaterial { - fn from(inner: Conductor) -> Self { +impl From> for GenericMaterial { + fn from(inner: Conductor) -> Self { Self::Conductor(inner) } } -impl From for GenericMaterial { - fn from(inner: LinearMagnet) -> Self { +impl From> for GenericMaterial { + fn from(inner: LinearMagnet) -> Self { Self::LinearMagnet(inner) } } -impl From for GenericMaterial { - fn from(inner: Pml) -> Self { +impl From> for GenericMaterial { + fn from(inner: Pml) -> Self { Self::Pml(inner) } } -impl From for GenericMaterial { - fn from(inner: Ferroxcube3R1) -> Self { +impl From> for GenericMaterial { + fn from(inner: Ferroxcube3R1) -> Self { Self::Ferroxcube3R1(inner) } } -impl From for GenericMaterial { - fn from(inner: MinimalSquare) -> Self { +impl From> for GenericMaterial { + fn from(inner: MinimalSquare) -> Self { Self::MinimalSquare(inner) } } -impl Material for GenericMaterial { - fn step_parameters_mut<'a>(&'a mut self) -> StepParametersMut<'a, flt::Real> { +impl Material for GenericMaterial { + fn step_parameters_mut<'a>(&'a mut self) -> StepParametersMut<'a, R> { use GenericMaterial::*; match self { Conductor(inner) => inner.step_parameters_mut(), @@ -445,7 +488,7 @@ impl Material for GenericMaterial { } } /// Return the magnetization. - fn m(&self) -> Vec3 { + fn m(&self) -> Vec3 { use GenericMaterial::*; match self { Conductor(inner) => inner.m(), @@ -456,7 +499,7 @@ impl Material for GenericMaterial { } } /// Called just before magnetic field is updated. Optionally change any internal state (e.g. magnetization). - fn step_b(&mut self, context: &CellState, delta_b: Vec3) { + fn step_b(&mut self, context: &CellState, delta_b: Vec3) { use GenericMaterial::*; match self { Conductor(inner) => inner.step_b(context, delta_b), @@ -470,27 +513,27 @@ impl Material for GenericMaterial { // #[enum_dispatch(Material)] #[derive(Clone, Serialize, Deserialize)] -pub enum GenericMaterialNoPml { - Conductor(Conductor), - LinearMagnet(LinearMagnet), - Ferroxcube3R1(Ferroxcube3R1), - MinimalSquare(MinimalSquare), +pub enum GenericMaterialNoPml { + Conductor(Conductor), + LinearMagnet(LinearMagnet), + Ferroxcube3R1(Ferroxcube3R1), + MinimalSquare(MinimalSquare), } -impl Default for GenericMaterialNoPml { +impl Default for GenericMaterialNoPml { fn default() -> Self { Conductor::default().into() } } -impl From for GenericMaterialNoPml { - fn from(inner: Conductor) -> Self { +impl From> for GenericMaterialNoPml { + fn from(inner: Conductor) -> Self { Self::Conductor(inner) } } -impl Material for GenericMaterialNoPml { - fn step_parameters_mut<'a>(&'a mut self) -> StepParametersMut<'a, flt::Real> { +impl Material for GenericMaterialNoPml { + fn step_parameters_mut<'a>(&'a mut self) -> StepParametersMut<'a, R> { use GenericMaterialNoPml::*; match self { Conductor(inner) => inner.step_parameters_mut(), @@ -500,7 +543,7 @@ impl Material for GenericMaterialNoPml { } } /// Return the magnetization. - fn m(&self) -> Vec3 { + fn m(&self) -> Vec3 { use GenericMaterialNoPml::*; match self { Conductor(inner) => inner.m(), @@ -510,7 +553,7 @@ impl Material for GenericMaterialNoPml { } } /// Called just before magnetic field is updated. Optionally change any internal state (e.g. magnetization). - fn step_b(&mut self, context: &CellState, delta_b: Vec3) { + fn step_b(&mut self, context: &CellState, delta_b: Vec3) { use GenericMaterialNoPml::*; match self { Conductor(inner) => inner.step_b(context, delta_b), @@ -525,26 +568,26 @@ impl Material for GenericMaterialNoPml { /// Materials which have only 1 Vec3. // #[enum_dispatch(Material)] #[derive(Clone, Serialize, Deserialize)] -pub enum GenericMaterialOneField { - Conductor(Conductor), - Ferroxcube3R1(Ferroxcube3R1), - MinimalSquare(MinimalSquare), +pub enum GenericMaterialOneField { + Conductor(Conductor), + Ferroxcube3R1(Ferroxcube3R1), + MinimalSquare(MinimalSquare), } -impl Default for GenericMaterialOneField { +impl Default for GenericMaterialOneField { fn default() -> Self { Conductor::default().into() } } -impl From for GenericMaterialOneField { - fn from(inner: Conductor) -> Self { +impl From> for GenericMaterialOneField { + fn from(inner: Conductor) -> Self { Self::Conductor(inner) } } -impl Material for GenericMaterialOneField { - fn step_parameters_mut<'a>(&'a mut self) -> StepParametersMut<'a, flt::Real> { +impl Material for GenericMaterialOneField { + fn step_parameters_mut<'a>(&'a mut self) -> StepParametersMut<'a, R> { use GenericMaterialOneField::*; match self { Conductor(inner) => inner.step_parameters_mut(), @@ -553,7 +596,7 @@ impl Material for GenericMaterialOneField { } } /// Return the magnetization. - fn m(&self) -> Vec3 { + fn m(&self) -> Vec3 { use GenericMaterialOneField::*; match self { Conductor(inner) => inner.m(), @@ -562,7 +605,7 @@ impl Material for GenericMaterialOneField { } } /// Called just before magnetic field is updated. Optionally change any internal state (e.g. magnetization). - fn step_b(&mut self, context: &CellState, delta_b: Vec3) { + fn step_b(&mut self, context: &CellState, delta_b: Vec3) { use GenericMaterialOneField::*; match self { Conductor(inner) => inner.step_b(context, delta_b), @@ -575,32 +618,32 @@ impl Material for GenericMaterialOneField { /// Database of common materials pub mod db { use super::*; - pub fn conductor(conductivity: R) -> Conductor { + pub fn conductor(conductivity: R2) -> Conductor { Conductor::new(conductivity) } - pub fn anisotropic_conductor(conductivity: Vec3) -> Conductor { + pub fn anisotropic_conductor(conductivity: Vec3) -> Conductor { Conductor::new_anisotropic(conductivity) } - pub fn copper() -> Conductor { + pub fn copper() -> Conductor { Conductor::new(50_000_000.0) } // See https://en.wikipedia.org/wiki/Permeability_(electromagnetism)#Values_for_some_common_materials /// This is a simplified form of iron annealed in H. - pub fn linear_annealed_iron() -> LinearMagnet { + pub fn linear_annealed_iron() -> LinearMagnet { LinearMagnet::new(200_000.0) } /// This is a simplified form of iron - pub fn linear_iron() -> LinearMagnet { + pub fn linear_iron() -> LinearMagnet { LinearMagnet::new(5000.0) } /// https://www.ferroxcube.com/upload/media/product/file/MDS/3r1.pdf - pub fn ferroxcube_3r1() -> Ferroxcube3R1 { + pub fn ferroxcube_3r1() -> Ferroxcube3R1 { Ferroxcube3R1::default() } - pub fn minimal_square_ferrite() -> MinimalSquare { + pub fn minimal_square_ferrite() -> MinimalSquare { MinimalSquare::default() } } @@ -608,7 +651,7 @@ pub mod db { #[cfg(test)] mod test { use super::*; - fn mh_curve_for_test() -> MHCurve { + fn mh_curve_for_test() -> MHCurve { MHCurve::new(&[ // rising Vec2::new( 10.0, 0.0), @@ -625,24 +668,21 @@ mod test { ]) } - fn assert_step_toward_symmetric(h: Flt, m: Flt, target_mh: Flt, target: Result<(Flt, Flt), (Flt, Flt)>) { + fn assert_step_toward_symmetric(h: f32, m: f32, target_mh: f32, target: Result, Vec2>) { let curve = mh_curve_for_test(); - let h = flt::Real::from_inner(h); - let m = flt::Real::from_inner(m); - let target_mh = flt::Real::from_inner(target_mh); let target = match target { - Ok((a, b)) => Ok((flt::Real::from_inner(a), flt::Real::from_inner(b))), - Err((a, b)) => Err((flt::Real::from_inner(a), flt::Real::from_inner(b))), + Ok(v) => Ok(v), + Err(v) => Err(v), }; let neg_target = match target { - Ok((a, b)) => Ok((-a, -b)), - Err((a, b)) => Err((-a, -b)), + Ok(v) => Ok(-v), + Err(v) => Err(-v), }; assert_eq!(curve.step_toward(h, m, target_mh), target); assert_eq!(curve.step_toward(-h, -m, -target_mh), neg_target); } - fn assert_move_to_symmetric(h: Flt, m: Flt, target_mh: Flt, target: (Flt, Flt)) { + fn assert_move_to_symmetric(h: f32, m: f32, target_mh: f32, target: (f32, f32)) { let curve = mh_curve_for_test(); assert_eq!(curve.move_to(h, m, target_mh), target); assert_eq!(curve.move_to(-h, -m, -target_mh), (-target.0, -target.1)); @@ -650,42 +690,42 @@ mod test { #[test] fn mh_curve_move_from_inner_to_inner() { - assert_step_toward_symmetric(0.0, 0.0, 5.0, Ok((5.0, 0.0))); - assert_step_toward_symmetric(0.0, 5.0, 10.0, Ok((5.0, 5.0))); + assert_step_toward_symmetric(0.0, 0.0, 5.0, Ok(Vec2::new(5.0, 0.0))); + assert_step_toward_symmetric(0.0, 5.0, 10.0, Ok(Vec2::new(5.0, 5.0))); - assert_step_toward_symmetric(-5.0, 5.0, -3.0, Ok((-8.0, 5.0))); - assert_step_toward_symmetric(-5.0, 5.0, 7.0, Ok((2.0, 5.0))); + assert_step_toward_symmetric(-5.0, 5.0, -3.0, Ok(Vec2::new(-8.0, 5.0))); + assert_step_toward_symmetric(-5.0, 5.0, 7.0, Ok(Vec2::new(2.0, 5.0))); - assert_step_toward_symmetric(5.0, -5.0, -3.0, Ok((2.0, -5.0))); - assert_step_toward_symmetric(5.0, -5.0, 3.0, Ok((8.0, -5.0))); + assert_step_toward_symmetric(5.0, -5.0, -3.0, Ok(Vec2::new(2.0, -5.0))); + assert_step_toward_symmetric(5.0, -5.0, 3.0, Ok(Vec2::new(8.0, -5.0))); } #[test] fn mh_curve_magnetize_along_edge() { // start of segment NOOP - assert_step_toward_symmetric(10.0, 0.0, 10.0, Ok((10.0, 0.0))); + assert_step_toward_symmetric(10.0, 0.0, 10.0, Ok(Vec2::new(10.0, 0.0))); // start of segment to middle of segment - assert_step_toward_symmetric(10.0, 0.0, 32.0, Ok((12.0, 20.0))); + assert_step_toward_symmetric(10.0, 0.0, 32.0, Ok(Vec2::new(12.0, 20.0))); // middle of segment NOOP - assert_step_toward_symmetric(12.0, 20.0, 32.0, Ok((12.0, 20.0))); + assert_step_toward_symmetric(12.0, 20.0, 32.0, Ok(Vec2::new(12.0, 20.0))); // middle of segment to middle of segment - assert_step_toward_symmetric(12.0, 20.0, 54.0, Ok((14.0, 40.0))); + assert_step_toward_symmetric(12.0, 20.0, 54.0, Ok(Vec2::new(14.0, 40.0))); // middle of segment to end of segment - assert_step_toward_symmetric(12.0, 20.0, 120.0, Err((20.0, 100.0))); + assert_step_toward_symmetric(12.0, 20.0, 120.0, Err(Vec2::new(20.0, 100.0))); } #[test] fn mh_curve_demagnetize_along_edge() { // start of segment NOOP - assert_step_toward_symmetric(30.0, 150.0, 180.0, Ok((30.0, 150.0))); + assert_step_toward_symmetric(30.0, 150.0, 180.0, Ok(Vec2::new(30.0, 150.0))); // start of segment to middle of segment - assert_step_toward_symmetric(30.0, 150.0, 160.0, Ok((20.0, 140.0))); + assert_step_toward_symmetric(30.0, 150.0, 160.0, Ok(Vec2::new(20.0, 140.0))); // middle of segment NOOP - assert_step_toward_symmetric(20.0, 140.0, 160.0, Ok((20.0, 140.0))); + assert_step_toward_symmetric(20.0, 140.0, 160.0, Ok(Vec2::new(20.0, 140.0))); // middle of segment to middle of segment - assert_step_toward_symmetric(20.0, 140.0, 140.0, Ok((10.0, 130.0))); + assert_step_toward_symmetric(20.0, 140.0, 140.0, Ok(Vec2::new(10.0, 130.0))); // middle of segment to end of segment - assert_step_toward_symmetric(20.0, 140.0, 120.0, Err((0.0, 120.0))); + assert_step_toward_symmetric(20.0, 140.0, 120.0, Err(Vec2::new(0.0, 120.0))); } #[test] diff --git a/src/post.rs b/src/post.rs index 2547671..4214807 100644 --- a/src/post.rs +++ b/src/post.rs @@ -111,7 +111,7 @@ 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>> = bincode::deserialize_from(reader)?; Ok(data.to_static()) })?; diff --git a/src/real.rs b/src/real.rs index bc4c361..59c2767 100644 --- a/src/real.rs +++ b/src/real.rs @@ -15,7 +15,7 @@ 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 + Copy + Clone + Default + Send + Sync { +pub trait Real: ToFloat + decorum_Real + IntrinsicOrd + AddAssign + MulAssign + fmt::LowerExp + fmt::Display + fmt::Debug + Copy + Clone + Default + Send + Sync + 'static { // TODO: fold with from_ fn from_primitive(p: P) -> Self { Self::from_f64(p.to_f64()) diff --git a/src/sim.rs b/src/sim.rs index 22c5856..ee3e9d7 100644 --- a/src/sim.rs +++ b/src/sim.rs @@ -1,7 +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, decorum_Real as _, Real as _, ToFloat as _, Zero as _}; +use crate::real::{self, Real as _, ToFloat as _}; use crate::stim::AbstractStimulus; use dyn_clone::{self, DynClone}; use log::trace; @@ -11,7 +11,7 @@ use serde::{Serialize, Deserialize}; use std::convert::From; use std::iter::Sum; -pub type StaticSim = SimState; +pub type StaticSim = SimState, Real>; #[derive(Default, Copy, Clone, PartialEq, Serialize, Deserialize)] pub struct PmlState { @@ -372,7 +372,7 @@ impl<'a> dyn GenericSim + 'a { } #[derive(Default, Clone, Serialize, Deserialize)] -pub struct SimState { +pub struct SimState, R=Real> { cells: Array3, e: Array3>, h: Array3>, @@ -1205,10 +1205,10 @@ mod test { let signal = [2.0, 0.0, 0.0]; // kernel: e(-t) // \int_0^1 e(-t) dt = [1 - exp(-1)] - let exp_neg_0 = (-0.0).exp(); - let exp_neg_1 = (-1.0).exp(); - let exp_neg_2 = (-2.0).exp(); - let exp_neg_3 = (-3.0).exp(); + let exp_neg_0 = (-0.0f64).exp(); + let exp_neg_1 = (-1.0f64).exp(); + let exp_neg_2 = (-2.0f64).exp(); + let exp_neg_3 = (-3.0f64).exp(); let expected = [ 2.0*(exp_neg_0 - exp_neg_1), 2.0*(exp_neg_1 - exp_neg_2), @@ -1223,10 +1223,10 @@ mod test { let signal = [2.0, 0.0, 0.0]; // kernel: e(-3*t) // \int_0^0.2 e(-3*t) dt = [1 - exp(-0.6)]/3 - let exp_neg_00 = (-0.0).exp(); - let exp_neg_06 = (-0.6).exp(); - let exp_neg_12 = (-1.2).exp(); - let exp_neg_18 = (-1.8).exp(); + let exp_neg_00 = (-0.0f64).exp(); + let exp_neg_06 = (-0.6f64).exp(); + let exp_neg_12 = (-1.2f64).exp(); + let exp_neg_18 = (-1.8f64).exp(); let expected = [ 2.0/3.0*(exp_neg_00 - exp_neg_06), 2.0/3.0*(exp_neg_06 - exp_neg_12), @@ -1241,10 +1241,10 @@ mod test { let signal = [2.0, 7.0, -3.0]; // kernel: e(-3*t) // \int_0^0.2 e(-3*t) dt = [1 - exp(-0.6)]/3 - let exp_neg_00 = (-0.0).exp(); - let exp_neg_06 = (-0.6).exp(); - let exp_neg_12 = (-1.2).exp(); - let exp_neg_18 = (-1.8).exp(); + let exp_neg_00 = (-0.0f64).exp(); + let exp_neg_06 = (-0.6f64).exp(); + let exp_neg_12 = (-1.2f64).exp(); + let exp_neg_18 = (-1.8f64).exp(); let expected = [ 2.0/3.0*(exp_neg_00 - exp_neg_06), 7.0/3.0*(exp_neg_00 - exp_neg_06) + 2.0/3.0*(exp_neg_06 - exp_neg_12), @@ -1261,10 +1261,10 @@ mod test { let signal = [2.0, 7.0, -3.0]; // kernel: 0.3 \delta(t) + 0.4 * e(-3*t) // \int_0^0.2 e(-3*t) dt = [1 - exp(-0.6)]/3 - let exp_neg_00 = (-0.0).exp(); - let exp_neg_06 = (-0.6).exp(); - let exp_neg_12 = (-1.2).exp(); - let exp_neg_18 = (-1.8).exp(); + let exp_neg_00 = (-0.0f64).exp(); + let exp_neg_06 = (-0.6f64).exp(); + let exp_neg_12 = (-1.2f64).exp(); + let exp_neg_18 = (-1.8f64).exp(); let expected_exp = [ 2.0/3.0*(exp_neg_00 - exp_neg_06), 7.0/3.0*(exp_neg_00 - exp_neg_06) + 2.0/3.0*(exp_neg_06 - exp_neg_12), @@ -1330,8 +1330,8 @@ mod test { // nabla_g = 0.1 df/dt + 2 f // Let f(t) = sin(2 t) // then f'(t) = 2 cos(2 t) - let f = Vec3::uniform((5.0).sin()); - let df_dt = Vec3::uniform(2.0 * (5.0).cos()); + let f = Vec3::uniform((5.0f64).sin()); + let df_dt = Vec3::uniform(2.0 * (5.0f64).cos()); let nabla_g = df_dt*0.1 + f*2.0; let actual_df_dt = solve_step_diff_eq( nabla_g, @@ -1353,9 +1353,9 @@ mod test { // Let f(t) = sin(2 t) // then f'(t) = 2 cos(2 t) // then \int f(t) = -1/2 cos(2 t) - let f = Vec3::uniform((5.0).sin()); - let df_dt = Vec3::uniform(2.0 * (5.0).cos()); - let f_int = Vec3::uniform(-0.5 * (5.0).cos()); + let f = Vec3::uniform((5.0f64).sin()); + let df_dt = Vec3::uniform(2.0 * (5.0f64).cos()); + let f_int = Vec3::uniform(-0.5 * (5.0f64).cos()); let nabla_g = df_dt*0.1 + f*2.0 + f_int*0.4; let actual_df_dt = solve_step_diff_eq( nabla_g, @@ -1378,10 +1378,10 @@ mod test { // then f'(t) = 2 cos(2 t) // then \int f(t) = -1/2 cos(2 t) // then \int \int f(t) = -1/4 sin(2 t) - let f = Vec3::unit()*(5.0).sin(); - let df_dt = Vec3::unit()*2.0 * (5.0).cos(); - let f_int = Vec3::unit()*-0.5 * (5.0).cos(); - let f_int_int = Vec3::unit()*-0.25 * (5.0).sin(); + let f = Vec3::unit()*(5.0f64).sin(); + let df_dt = Vec3::unit()*2.0 * (5.0f64).cos(); + let f_int = Vec3::unit()*-0.5 * (5.0f64).cos(); + let f_int_int = Vec3::unit()*-0.25 * (5.0f64).sin(); let nabla_g = df_dt*0.1 + f*2.0 + f_int*0.4 + f_int_int*0.3; let actual_df_dt = solve_step_diff_eq( nabla_g, @@ -1569,7 +1569,7 @@ 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) { + 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 { @@ -1616,7 +1616,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> { let mut state = SimState::new(size, 1e-6); let timestep = state.timestep(); state.fill_boundary_using(size/4, |boundary_ness| { @@ -1782,7 +1782,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> { let mut state = SimState::new(size, 1e-6); let timestep = state.timestep(); state.fill_boundary_using(size/4, |boundary_ness| { @@ -1805,7 +1805,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);