Parameterize most of SimState and Material over Real.

Some things still work only for coremem::Real. Need to troubleshoot
those.

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::<GenericMaterial>::new(Index::new(size, size, size), 1e-5);
+            let mut driver = Driver::<GenericMaterial<Real>>::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::<GenericMaterialNoPml>::new(Index::new(size, size, size), 1e-5);
+            let mut driver = Driver::<GenericMaterialNoPml<Real>>::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::<GenericMaterialOneField>::new(Index::new(size, size, size), 1e-5);
+            let mut driver = Driver::<GenericMaterialOneField<Real>>::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::<GenericMaterial>::new(Index::new(size, size, size), 1e-5);
+            let mut driver = Driver::<GenericMaterial<Real>>::new(Index::new(size, size, size), 1e-5);
             driver.add_pml_boundary(Index::new(thickness, thickness, thickness));
             b.iter(|| driver.step())
         });

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<mat::GenericMaterial> = Driver::new(size_px, feat_size);
+    let mut driver: Driver<mat::GenericMaterial<Real>> = Driver::new(size_px, feat_size);
     driver.set_steps_per_frame(400);
     driver.set_steps_per_stim(1);
     let base = "minimal_torus-6";

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<mat::GenericMaterial> = Driver::new(size_px, feat_size);
+        let mut driver: Driver<mat::GenericMaterial<Real>> = Driver::new(size_px, feat_size);
         //driver.set_steps_per_frame(8);
         //driver.set_steps_per_frame(20);
         //driver.set_steps_per_frame(40);

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<mat::GenericMaterial> = Driver::new(size_px, feat_size);
+    let mut driver: Driver<mat::GenericMaterial<Real>> = 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";

src/bin/diagnostics.rs

@@ -1,7 +1,7 @@
 use coremem::mat::*;
 
 fn main() {
-    let m3r1 = Ferroxcube3R1::curve();
+    let m3r1 = Ferroxcube3R1::<f32>::curve();
     let extremes = m3r1.extremes();
     println!("3R1 Extremes: H={}, M={}", extremes.x(), extremes.y());
 }

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<M=GenericMaterial, R=flt::Real> {
+pub struct Driver<M=GenericMaterial<flt::Real>, R=flt::Real> {
     pub state: SimState<M, R>,
     renderer: Arc<MultiRenderer<SimState<M, R>>>,
     // TODO: use Rayon's thread pool?
@@ -192,7 +192,7 @@ impl<R: Real + Send + Sync + 'static, M: Material<R> + Clone + Default + Send +
     }
 }
 
-impl<R: Real, M: Material<R> + From<mat::Pml>> Driver<M, R> {
+impl<R: Real, M: Material<R> + From<mat::Pml<R>>> Driver<M, R> {
     pub fn add_pml_boundary<C: Coord>(&mut self, thickness: C) {
         let timestep = self.state.timestep();
         self.state.fill_boundary_using(thickness, |boundary_ness| {
@@ -203,7 +203,7 @@ impl<R: Real, M: Material<R> + From<mat::Pml>> Driver<M, R> {
     }
 }
 
-impl<R: Real, M: Material<R> + From<mat::Conductor>> Driver<M, R> {
+impl<R: Real, M: Material<R> + From<mat::Conductor<R>>> Driver<M, R> {
     pub fn add_classical_boundary<C: Coord>(&mut self, thickness: C) {
         let timestep = self.state.timestep();
         self.state.fill_boundary_using(thickness, |boundary_ness| {

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<R=f32> {
     pub x: R,
     pub y: R,

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<R: Real = flt::Real> {
+pub trait Material<R: Real> {
     fn step_parameters_mut<'a>(&'a mut self) -> StepParametersMut<'a, R> {
         // by default, behave as a vacuum
         StepParametersMut::default()
@@ -42,14 +45,14 @@ impl<R: Real, M: Material<R>> MaterialExt<R> 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 struct Static<R> {
-    pub conductivity: Vec3<flt::Real>,
+    pub conductivity: Vec3<R>,
     // pub pml: Option<(PmlState, PmlParameters)>,
-    pub m: Vec3<flt::Real>,
+    pub m: Vec3<R>,
 }
 
-impl Static {
+impl<R: Real> Static<R> {
-    pub fn from_material<M: Material<flt::Real>>(m: &M) -> Self {
+    pub fn from_material<M: Material<R>>(m: &M) -> Self {
         let p = m.step_parameters();
         Self {
             conductivity: p.conductivity(),
@@ -63,20 +66,20 @@ impl Static {
     // }
 }
 
-impl Material<flt::Real> for Static {
+impl<R: Real> Material<R> for Static<R> {
-    fn step_parameters_mut<'a>(&'a mut self) -> StepParametersMut<'a, flt::Real> {
+    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<flt::Real> {
+    fn m(&self) -> Vec3<R> {
         self.m
     }
 }
 
-impl<T> From<T> for Static
+impl<R: Real, T> From<T> for Static<R>
-where T: Into<GenericMaterial>
+where T: Into<GenericMaterial<R>>
 {
     fn from(mat: T) -> Self {
         let generic = mat.into();
@@ -86,45 +89,45 @@ where T: Into<GenericMaterial>
 
 /// Material which has a conductivity parameter, but cannot be magnetized
 #[derive(Clone, Default, Serialize, Deserialize)]
-pub struct Conductor {
+pub struct Conductor<R> {
-    conductivity: Vec3<flt::Real>,
+    conductivity: Vec3<R>,
 }
 
-impl Conductor {
+impl<R: Real> Conductor<R> {
-    pub fn new<R: Real>(conductivity: R) -> Self {
+    pub fn new<R2: Real>(conductivity: R2) -> Self {
         Self {
             conductivity: Vec3::uniform(conductivity).cast()
         }
     }
-    pub fn new_anisotropic<R: Real>(conductivity: Vec3<R>) -> Self {
+    pub fn new_anisotropic<R2: Real>(conductivity: Vec3<R2>) -> Self {
         Self {
             conductivity: conductivity.cast(),
         }
     }
 }
 
-impl Material<flt::Real> for Conductor {
+impl<R: Real> Material<R> for Conductor<R> {
-    fn step_parameters_mut<'a>(&'a mut self) -> StepParametersMut<'a, flt::Real> {
+    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<R> {
     /// \mu_r
-    relative_permeability: Vec3<flt::Real>,
+    relative_permeability: Vec3<R>,
-    m: Vec3<flt::Real>,
+    m: Vec3<R>,
 }
 
-impl LinearMagnet {
+impl<R: Real> LinearMagnet<R> {
-    pub fn new<R: Real>(relative_permeability: R) -> Self {
+    pub fn new<R2: Real>(relative_permeability: R2) -> Self {
         Self {
             relative_permeability: Vec3::uniform(relative_permeability).cast(),
             m: Vec3::zero(),
         }
     }
-    pub fn new_anisotropic<R: Real>(relative_permeability: Vec3<R>) -> Self {
+    pub fn new_anisotropic<R2: Real>(relative_permeability: Vec3<R2>) -> Self {
         Self {
             relative_permeability: relative_permeability.cast(),
             m: Vec3::zero()
@@ -132,11 +135,11 @@ impl LinearMagnet {
     }
 }
 
-impl Material<flt::Real> for LinearMagnet {
+impl<R: Real> Material<R> for LinearMagnet<R> {
-    fn m(&self) -> Vec3<flt::Real> {
+    fn m(&self) -> Vec3<R> {
         self.m
     }
-    fn step_b(&mut self, _context: &CellState<flt::Real>, delta_b: Vec3<flt::Real>) {
+    fn step_b(&mut self, _context: &CellState<R>, delta_b: Vec3<R>) {
         //```tex
         // $B = \mu_0 (H + M) = \mu_0 \mu_r H$
         // $\mu_r H = H + M$
@@ -145,80 +148,107 @@ impl Material<flt::Real> 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<flt::Real>, PmlParameters<flt::Real>);
+pub struct Pml<R>(PmlState<R>, PmlParameters<R>);
 
-impl Pml {
+impl<R: Real> Pml<R> {
-    pub fn new<R: Real>(pseudo_conductivity: Vec3<R>) -> Self {
+    pub fn new<R2: Real>(pseudo_conductivity: Vec3<R2>) -> Self {
         Self(PmlState::new(), PmlParameters::new(pseudo_conductivity))
     }
 }
 
-impl Material<flt::Real> for Pml {
+impl<R: Real> Material<R> for Pml<R> {
-    fn step_parameters_mut<'a>(&'a mut self) -> StepParametersMut<'a, flt::Real> {
+    fn step_parameters_mut<'a>(&'a mut self) -> StepParametersMut<'a, R> {
         StepParametersMut::default().with_pml(&mut self.0, self.1)
     }
 }
 
-pub trait PiecewiseLinearFerromagnet {
+// pub trait PiecewiseLinearFerromagnet<R> {
-    fn curve() -> &'static MHCurve;
+//     fn curve() -> &'static MHCurve<R>;
-    fn conductivity() -> Flt;
+//     fn conductivity() -> R;
-    fn m(&self) -> Vec3<flt::Real>;
+//     fn m(&self) -> Vec3<R>;
-    fn m_mut(&mut self) -> &mut Vec3<flt::Real>;
+//     fn m_mut(&mut self) -> &mut Vec3<R>;
-}
+// }
+// 
+// impl<R: Real, T: PiecewiseLinearFerromagnet<R>> Material<R> 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<R> {
+//         self.m()
+//     }
+//     fn step_b(&mut self, context: &CellState<R>, delta_b: Vec3<R>) {
+//         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<T: PiecewiseLinearFerromagnet> Material<flt::Real> for T {
+fn step_linear_ferro<R: Real>(m_mut: &mut Vec3<R>, mh_curve: &MHCurve<R>, context: &CellState<R>, delta_b: Vec3<R>) {
-    fn step_parameters_mut<'a>(&'a mut self) -> StepParametersMut<'a, flt::Real> {
+    trace!("step_b enter");
-        let c = T::conductivity();
+    let (h, m) = (context.h(), *m_mut);
-        StepParametersMut::default().with_conductivity(Vec3::uniform(c))
+    let target_hm = h + m + delta_b * R::mu0_inv();
-    }
-    fn m(&self) -> Vec3<flt::Real> {
-        self.m()
-    }
-    fn step_b(&mut self, context: &CellState<flt::Real>, delta_b: Vec3<flt::Real>) {
-        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();
 
-        // TODO: this is probably not the best way to generalize a BH curve into 3d.
+    // TODO: this is probably not the best way to generalize a BH curve into 3d.
-        let (_hx, mx) = mh_curve.move_to(
+    let (_hx, mx) = mh_curve.move_to(
-            Flt::from_primitive(h.x()),
+        h.x(),
-            Flt::from_primitive(m.x()),
+        m.x(),
-            Flt::from_primitive(target_hm.x())
+        target_hm.x(),
-        );
+    );
-        let (_hy, my) = mh_curve.move_to(
+    let (_hy, my) = mh_curve.move_to(
-            Flt::from_primitive(h.y()),
+        h.y(),
-            Flt::from_primitive(m.y()),
+        m.y(),
-            Flt::from_primitive(target_hm.y())
+        target_hm.y(),
-        );
+    );
-        let (_hz, mz) = mh_curve.move_to(
+    let (_hz, mz) = mh_curve.move_to(
-            Flt::from_primitive(h.z()),
+        h.z(),
-            Flt::from_primitive(m.z()),
+        m.z(),
-            Flt::from_primitive(target_hm.z())
+        target_hm.z(),
-        );
+    );
 
-        *self.m_mut() = Vec3::new(mx, my, mz).cast();
+    *m_mut = Vec3::new(mx, my, mz);
-        // let ret = Vec3::new(hx, hy, hz);
+    // let ret = Vec3::new(hx, hy, hz);
-        trace!("step_b end");
+    trace!("step_b end");
-    }
 }
 
 /// M as a function of H
 #[derive(Clone)]
-pub struct MHCurve {
+pub struct MHCurve<R> {
-    geom: Polygon2d<flt::Real>,
+    geom: Polygon2d<R>,
 }
 
-impl MHCurve {
+impl<R: Real> MHCurve<R> {
     /// 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<R: Real>(points: &[Vec2<R>]) -> Self {
+    fn new<R2: Real>(points: &[Vec2<R2>]) -> Self {
         let full_pts: Vec<_> =
             points.iter().cloned()
             .chain(points.iter().cloned().map(|p| -p))
@@ -230,15 +260,15 @@ impl MHCurve {
         }
     }
 
-    fn from_bh<R: Real>(points: &[(R, R)]) -> Self {
+    fn from_bh<R2: Real>(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<R: Real>(points: &[(R, R)]) -> Self {
+    fn from_mh<R2: Real>(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<flt::Real> {
+    pub fn extremes(&self) -> Vec2<R> {
         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<R>, Vec2<R>> {
         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::<Vec<_>>());
             });
             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()),
+                Ok(v) => break (v.x(), v.y()),
-                Err((x, y)) => {
+                Err(v) => {
-                    h = x;
+                    h = v.x();
-                    m = y;
+                    m = v.y();
                 },
             }
             if i % 2048 == 0 {
@@ -324,14 +353,18 @@ impl MHCurve {
 }
 
 #[derive(Default, Copy, Clone, Serialize, Deserialize)]
-pub struct Ferroxcube3R1 {
+pub struct Ferroxcube3R1<R> {
-    m: Vec3<flt::Real>,
+    m: Vec3<R>,
 }
 
-impl PiecewiseLinearFerromagnet for Ferroxcube3R1 {
+impl<R: Real> Ferroxcube3R1<R> {
-    fn curve() -> &'static MHCurve {
+    pub fn curve() -> &'static MHCurve<R> {
         lazy_static! {
-            static ref FERROXCUBE_3R1: MHCurve = MHCurve::from_bh(&[
+            static ref curves: Mutex<HashMap<TypeId, Box<dyn Any + Send>>> = Mutex::new(HashMap::new());
+        }
+        let mut lock = curves.lock().unwrap();
+        let curve = lock.entry(TypeId::of::<R>()).or_insert_with(|| {
+            Box::new(MHCurve::<R>::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),
-            ]);
+            ]))
-        }
+        }).downcast_ref::<MHCurve<R>>().unwrap();
-        &*FERROXCUBE_3R1
+        unsafe { std::mem::transmute::<&MHCurve<R>, &'static MHCurve<R>>(curve) }
     }
-    fn conductivity() -> Flt {
+}
-        1e-3
+
+impl<R: Real> Material<R> for Ferroxcube3R1<R> {
+    fn step_b(&mut self, context: &CellState<R>, delta_b: Vec3<R>) {
+        step_linear_ferro(&mut self.m, Self::curve(), context, delta_b)
     }
-    fn m(&self) -> Vec3<flt::Real> {
+    fn m(&self) -> Vec3<R> {
         self.m
     }
-    fn m_mut(&mut self) -> &mut Vec3<flt::Real> {
+    fn step_parameters_mut<'a>(&'a mut self) -> StepParametersMut<'a, R> {
-        &mut self.m
+        StepParametersMut::default().with_conductivity(Vec3::uniform(1e-3))
     }
 }
 
 /// Simple, square-loop ferrite
 #[derive(Default, Copy, Clone, Serialize, Deserialize)]
-pub struct MinimalSquare {
+pub struct MinimalSquare<R> {
-    m: Vec3<flt::Real>,
+    m: Vec3<R>,
 }
 
-impl PiecewiseLinearFerromagnet for MinimalSquare {
+impl<R: Real> MinimalSquare<R> {
-    fn curve() -> &'static MHCurve {
+    pub fn curve() -> &'static MHCurve<R> {
         lazy_static! {
-            static ref CURVE: MHCurve = MHCurve::from_mh(&[
+            static ref curves: Mutex<HashMap<TypeId, Box<dyn Any + Send>>> = Mutex::new(HashMap::new());
+        }
+        let mut lock = curves.lock().unwrap();
+        let curve = lock.entry(TypeId::of::<R>()).or_insert_with(|| {
+            Box::new(MHCurve::<R>::from_bh(&[
                 (  1.0,       0.0),
                 (  2.0, 1000000.0),
                 // Falling
                 (  0.0,  900000.0),
-            ]);
+            ]))
-        }
+        }).downcast_ref::<MHCurve<R>>().unwrap();
-        &*CURVE
+        unsafe { std::mem::transmute::<&MHCurve<R>, &'static MHCurve<R>>(curve) }
     }
-    fn conductivity() -> Flt {
+}
-        1e-3
+
+impl<R: Real> Material<R> for MinimalSquare<R> {
+    fn step_b(&mut self, context: &CellState<R>, delta_b: Vec3<R>) {
+        step_linear_ferro(&mut self.m, Self::curve(), context, delta_b)
     }
-    fn m(&self) -> Vec3<flt::Real> {
+    fn m(&self) -> Vec3<R> {
         self.m
     }
-    fn m_mut(&mut self) -> &mut Vec3<flt::Real> {
+    fn step_parameters_mut<'a>(&'a mut self) -> StepParametersMut<'a, R> {
-        &mut self.m
+        StepParametersMut::default().with_conductivity(Vec3::uniform(1e-3))
     }
 }
 
 
 // #[enum_dispatch(Material)]
 #[derive(Clone, Serialize, Deserialize)]
-pub enum GenericMaterial {
+pub enum GenericMaterial<R> {
-    Conductor(Conductor),
+    Conductor(Conductor<R>),
-    LinearMagnet(LinearMagnet),
+    LinearMagnet(LinearMagnet<R>),
-    Pml(Pml),
+    Pml(Pml<R>),
-    Ferroxcube3R1(Ferroxcube3R1),
+    Ferroxcube3R1(Ferroxcube3R1<R>),
-    MinimalSquare(MinimalSquare),
+    MinimalSquare(MinimalSquare<R>),
 }
 
-impl Default for GenericMaterial {
+impl<R: Real> Default for GenericMaterial<R> {
     fn default() -> Self {
         Conductor::default().into()
     }
 }
 
-impl From<Conductor> for GenericMaterial {
+impl<R> From<Conductor<R>> for GenericMaterial<R> {
-    fn from(inner: Conductor) -> Self {
+    fn from(inner: Conductor<R>) -> Self {
         Self::Conductor(inner)
     }
 }
 
-impl From<LinearMagnet> for GenericMaterial {
+impl<R> From<LinearMagnet<R>> for GenericMaterial<R> {
-    fn from(inner: LinearMagnet) -> Self {
+    fn from(inner: LinearMagnet<R>) -> Self {
         Self::LinearMagnet(inner)
     }
 }
 
-impl From<Pml> for GenericMaterial {
+impl<R> From<Pml<R>> for GenericMaterial<R> {
-    fn from(inner: Pml) -> Self {
+    fn from(inner: Pml<R>) -> Self {
         Self::Pml(inner)
     }
 }
 
-impl From<Ferroxcube3R1> for GenericMaterial {
+impl<R> From<Ferroxcube3R1<R>> for GenericMaterial<R> {
-    fn from(inner: Ferroxcube3R1) -> Self {
+    fn from(inner: Ferroxcube3R1<R>) -> Self {
         Self::Ferroxcube3R1(inner)
     }
 }
 
-impl From<MinimalSquare> for GenericMaterial {
+impl<R> From<MinimalSquare<R>> for GenericMaterial<R> {
-    fn from(inner: MinimalSquare) -> Self {
+    fn from(inner: MinimalSquare<R>) -> Self {
         Self::MinimalSquare(inner)
     }
 }
 
-impl Material<flt::Real> for GenericMaterial {
+impl<R: Real> Material<R> for GenericMaterial<R> {
-    fn step_parameters_mut<'a>(&'a mut self) -> StepParametersMut<'a, flt::Real> {
+    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<flt::Real> for GenericMaterial {
         }
     }
     /// Return the magnetization.
-    fn m(&self) -> Vec3<flt::Real> {
+    fn m(&self) -> Vec3<R> {
         use GenericMaterial::*;
         match self {
             Conductor(inner) => inner.m(),
@@ -456,7 +499,7 @@ impl Material<flt::Real> for GenericMaterial {
         }
     }
     /// Called just before magnetic field is updated. Optionally change any internal state (e.g. magnetization).
-    fn step_b(&mut self, context: &CellState<flt::Real>, delta_b: Vec3<flt::Real>) {
+    fn step_b(&mut self, context: &CellState<R>, delta_b: Vec3<R>) {
         use GenericMaterial::*;
         match self {
             Conductor(inner) => inner.step_b(context, delta_b),
@@ -470,27 +513,27 @@ impl Material<flt::Real> for GenericMaterial {
 
 // #[enum_dispatch(Material)]
 #[derive(Clone, Serialize, Deserialize)]
-pub enum GenericMaterialNoPml {
+pub enum GenericMaterialNoPml<R> {
-    Conductor(Conductor),
+    Conductor(Conductor<R>),
-    LinearMagnet(LinearMagnet),
+    LinearMagnet(LinearMagnet<R>),
-    Ferroxcube3R1(Ferroxcube3R1),
+    Ferroxcube3R1(Ferroxcube3R1<R>),
-    MinimalSquare(MinimalSquare),
+    MinimalSquare(MinimalSquare<R>),
 }
 
-impl Default for GenericMaterialNoPml {
+impl<R: Real> Default for GenericMaterialNoPml<R> {
     fn default() -> Self {
         Conductor::default().into()
     }
 }
 
-impl From<Conductor> for GenericMaterialNoPml {
+impl<R> From<Conductor<R>> for GenericMaterialNoPml<R> {
-    fn from(inner: Conductor) -> Self {
+    fn from(inner: Conductor<R>) -> Self {
         Self::Conductor(inner)
     }
 }
 
-impl Material<flt::Real> for GenericMaterialNoPml {
+impl<R: Real> Material<R> for GenericMaterialNoPml<R> {
-    fn step_parameters_mut<'a>(&'a mut self) -> StepParametersMut<'a, flt::Real> {
+    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<flt::Real> for GenericMaterialNoPml {
         }
     }
     /// Return the magnetization.
-    fn m(&self) -> Vec3<flt::Real> {
+    fn m(&self) -> Vec3<R> {
         use GenericMaterialNoPml::*;
         match self {
             Conductor(inner) => inner.m(),
@@ -510,7 +553,7 @@ impl Material<flt::Real> for GenericMaterialNoPml {
         }
     }
     /// Called just before magnetic field is updated. Optionally change any internal state (e.g. magnetization).
-    fn step_b(&mut self, context: &CellState<flt::Real>, delta_b: Vec3<flt::Real>) {
+    fn step_b(&mut self, context: &CellState<R>, delta_b: Vec3<R>) {
         use GenericMaterialNoPml::*;
         match self {
             Conductor(inner) => inner.step_b(context, delta_b),
@@ -525,26 +568,26 @@ impl Material<flt::Real> for GenericMaterialNoPml {
 /// Materials which have only 1 Vec3.
 // #[enum_dispatch(Material)]
 #[derive(Clone, Serialize, Deserialize)]
-pub enum GenericMaterialOneField {
+pub enum GenericMaterialOneField<R> {
-    Conductor(Conductor),
+    Conductor(Conductor<R>),
-    Ferroxcube3R1(Ferroxcube3R1),
+    Ferroxcube3R1(Ferroxcube3R1<R>),
-    MinimalSquare(MinimalSquare),
+    MinimalSquare(MinimalSquare<R>),
 }
 
-impl Default for GenericMaterialOneField {
+impl<R: Real> Default for GenericMaterialOneField<R> {
     fn default() -> Self {
         Conductor::default().into()
     }
 }
 
-impl From<Conductor> for GenericMaterialOneField {
+impl<R> From<Conductor<R>> for GenericMaterialOneField<R> {
-    fn from(inner: Conductor) -> Self {
+    fn from(inner: Conductor<R>) -> Self {
         Self::Conductor(inner)
     }
 }
 
-impl Material<flt::Real> for GenericMaterialOneField {
+impl<R: Real> Material<R> for GenericMaterialOneField<R> {
-    fn step_parameters_mut<'a>(&'a mut self) -> StepParametersMut<'a, flt::Real> {
+    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<flt::Real> for GenericMaterialOneField {
         }
     }
     /// Return the magnetization.
-    fn m(&self) -> Vec3<flt::Real> {
+    fn m(&self) -> Vec3<R> {
         use GenericMaterialOneField::*;
         match self {
             Conductor(inner) => inner.m(),
@@ -562,7 +605,7 @@ impl Material<flt::Real> for GenericMaterialOneField {
         }
     }
     /// Called just before magnetic field is updated. Optionally change any internal state (e.g. magnetization).
-    fn step_b(&mut self, context: &CellState<flt::Real>, delta_b: Vec3<flt::Real>) {
+    fn step_b(&mut self, context: &CellState<R>, delta_b: Vec3<R>) {
         use GenericMaterialOneField::*;
         match self {
             Conductor(inner) => inner.step_b(context, delta_b),
@@ -575,32 +618,32 @@ impl Material<flt::Real> for GenericMaterialOneField {
 /// Database of common materials
 pub mod db {
     use super::*;
-    pub fn conductor<R: Real>(conductivity: R) -> Conductor {
+    pub fn conductor<R: Real, R2: Real>(conductivity: R2) -> Conductor<R> {
         Conductor::new(conductivity)
     }
-    pub fn anisotropic_conductor<R: Real>(conductivity: Vec3<R>) -> Conductor {
+    pub fn anisotropic_conductor<R: Real, R2: Real>(conductivity: Vec3<R2>) -> Conductor<R> {
         Conductor::new_anisotropic(conductivity)
     }
 
-    pub fn copper() -> Conductor {
+    pub fn copper<R: Real>() -> Conductor<R> {
         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<R: Real>() -> LinearMagnet<R> {
         LinearMagnet::new(200_000.0)
     }
     /// This is a simplified form of iron
-    pub fn linear_iron() -> LinearMagnet {
+    pub fn linear_iron<R: Real>() -> LinearMagnet<R> {
         LinearMagnet::new(5000.0)
     }
 
     /// https://www.ferroxcube.com/upload/media/product/file/MDS/3r1.pdf
-    pub fn ferroxcube_3r1() -> Ferroxcube3R1 {
+    pub fn ferroxcube_3r1<R: Real>() -> Ferroxcube3R1<R> {
         Ferroxcube3R1::default()
     }
-    pub fn minimal_square_ferrite() -> MinimalSquare {
+    pub fn minimal_square_ferrite<R: Real>() -> MinimalSquare<R> {
         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<f32> {
         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<f32>, Vec2<f32>>) {
         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))),
+            Ok(v) => Ok(v),
-            Err((a, b)) => Err((flt::Real::from_inner(a), flt::Real::from_inner(b))),
+            Err(v) => Err(v),
         };
         let neg_target = match target {
-            Ok((a, b)) => Ok((-a, -b)),
+            Ok(v) => Ok(-v),
-            Err((a, b)) => Err((-a, -b)),
+            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, 0.0, 5.0, Ok(Vec2::new(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, 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, -3.0, Ok(Vec2::new(-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, 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(Vec2::new(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(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]

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<SimState<GenericMaterial>> =
+            let data: SerializedFrame<SimState<GenericMaterial<crate::flt::Real>>> =
                 bincode::deserialize_from(reader)?;
             Ok(data.to_static())
         })?;

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_<blah>
     fn from_primitive<P: ToFloat>(p: P) -> Self {
         Self::from_f64(p.to_f64())

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<mat::Static, Real>;
+pub type StaticSim = SimState<mat::Static<Real>, Real>;
 
 #[derive(Default, Copy, Clone, PartialEq, Serialize, Deserialize)]
 pub struct PmlState<R> {
@@ -372,7 +372,7 @@ impl<'a> dyn GenericSim + 'a {
 }
 
 #[derive(Default, Clone, Serialize, Deserialize)]
-pub struct SimState<M=GenericMaterial, R=Real> {
+pub struct SimState<M=GenericMaterial<Real>, R=Real> {
     cells: Array3<M>,
     e: Array3<Vec3<R>>,
     h: Array3<Vec3<R>>,
@@ -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_0 = (-0.0f64).exp();
-        let exp_neg_1 = (-1.0).exp();
+        let exp_neg_1 = (-1.0f64).exp();
-        let exp_neg_2 = (-2.0).exp();
+        let exp_neg_2 = (-2.0f64).exp();
-        let exp_neg_3 = (-3.0).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_00 = (-0.0f64).exp();
-        let exp_neg_06 = (-0.6).exp();
+        let exp_neg_06 = (-0.6f64).exp();
-        let exp_neg_12 = (-1.2).exp();
+        let exp_neg_12 = (-1.2f64).exp();
-        let exp_neg_18 = (-1.8).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_00 = (-0.0f64).exp();
-        let exp_neg_06 = (-0.6).exp();
+        let exp_neg_06 = (-0.6f64).exp();
-        let exp_neg_12 = (-1.2).exp();
+        let exp_neg_12 = (-1.2f64).exp();
-        let exp_neg_18 = (-1.8).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_00 = (-0.0f64).exp();
-        let exp_neg_06 = (-0.6).exp();
+        let exp_neg_06 = (-0.6f64).exp();
-        let exp_neg_12 = (-1.2).exp();
+        let exp_neg_12 = (-1.2f64).exp();
-        let exp_neg_18 = (-1.8).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 f = Vec3::uniform((5.0f64).sin());
-        let df_dt = Vec3::uniform(2.0 * (5.0).cos());
+        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 f = Vec3::uniform((5.0f64).sin());
-        let df_dt = Vec3::uniform(2.0 * (5.0).cos());
+        let df_dt = Vec3::uniform(2.0 * (5.0f64).cos());
-        let f_int = Vec3::uniform(-0.5 * (5.0).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 f = Vec3::unit()*(5.0f64).sin();
-        let df_dt = Vec3::unit()*2.0 * (5.0).cos();
+        let df_dt = Vec3::unit()*2.0 * (5.0f64).cos();
-        let f_int = Vec3::unit()*-0.5 * (5.0).cos();
+        let f_int = Vec3::unit()*-0.5 * (5.0f64).cos();
-        let f_int_int = Vec3::unit()*-0.25 * (5.0).sin();
+        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<M: Into<mat::Static>>(mat: M) -> (f32, f32) {
+    fn conductor_test<M: Into<mat::Static<Real>>>(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<Pml> {
+    fn state_for_pml(size: Index) -> SimState<Pml<Real>> {
         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<Pml> {
+    fn state_for_monodirectional_pml(size: Index) -> SimState<Pml<Real>> {
         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<F: Fn(Vec3<f32>) -> Vec3<f32>>(
         size: Index, e: Vec3<f32>, shuffle: F
     ) {
-        let mut state = SimState::new(size, 1e-6);
+        let mut state = SimState::<Pml<Real>>::new(size, 1e-6);
         let timestep = state.timestep();
         state.fill_boundary_using(size/4, |boundary_ness| {
             let b = boundary_ness.elem_pow(3.0);