app: multi_core_inverter: more precise clock management

try to control the edges when the clock is release to prevent ringing.

Cargo.lock

@@ -160,6 +160,7 @@ dependencies = [
 name = "buffer_proto5"
 version = "0.1.0"
 dependencies = [
+ "bincode",
  "coremem",
  "log",
  "serde",
@@ -311,13 +312,12 @@ version = "0.1.0"
 dependencies = [
  "bincode",
  "common_macros",
- "coremem_types",
+ "coremem_cross",
  "criterion",
  "crossterm",
  "csv",
  "dashmap",
  "dyn-clone",
- "enum_dispatch",
  "env_logger",
  "float_eq",
  "font8x8",
@@ -344,6 +344,13 @@ dependencies = [
  "y4m",
 ]
 
+[[package]]
+name = "coremem_cross"
+version = "0.2.0"
+dependencies = [
+ "serde",
+]
+
 [[package]]
 name = "coremem_post"
 version = "0.1.0"
@@ -358,13 +365,6 @@ dependencies = [
  "structopt",
 ]
 
-[[package]]
-name = "coremem_types"
-version = "0.1.0"
-dependencies = [
- "serde",
-]
-
 [[package]]
 name = "crc32fast"
 version = "1.3.2"
@@ -577,18 +577,6 @@ version = "1.7.0"
 source = "registry+https://github.com/rust-lang/crates.io-index"
 checksum = "3f107b87b6afc2a64fd13cac55fe06d6c8859f12d4b14cbcdd2c67d0976781be"
 
-[[package]]
-name = "enum_dispatch"
-version = "0.3.8"
-source = "registry+https://github.com/rust-lang/crates.io-index"
-checksum = "0eb359f1476bf611266ac1f5355bc14aeca37b299d0ebccc038ee7058891c9cb"
-dependencies = [
- "once_cell",
- "proc-macro2",
- "quote",
- "syn",
-]
-
 [[package]]
 name = "env_logger"
 version = "0.9.0"
@@ -2095,7 +2083,7 @@ source = "git+https://github.com/EmbarkStudios/rust-gpu#0866cf591a7fdbbd15bdb346
 name = "spirv_backend"
 version = "0.1.0"
 dependencies = [
- "coremem_types",
+ "coremem_cross",
  "spirv-std",
 ]
 

Cargo.toml

@@ -5,7 +5,7 @@ members = [
     "crates/spirv_backend",
     "crates/spirv_backend_builder",
     "crates/spirv_backend_runner",
-    "crates/types",
+    "crates/cross",
     "crates/post",
 
     "crates/applications/buffer_proto5",

crates/applications/archive/buffer_proto1.rs

@@ -60,12 +60,12 @@ fn main() {
     driver.set_steps_per_stim(1000);
     //driver.fill_region(&ferro1_region, mat::db::linear_iron());
     // Original, 3R1-LIKE ferromagnet (only a vague likeness), sr-latch-8:
-    // driver.fill_region(&ferro1_region, mat::MBFerromagnet::new(-0.3899, 0.3900, 310_000.0));
-    // driver.fill_region(&ferro2_region, mat::MBFerromagnet::new(-0.3899, 0.3900, 310_000.0));
+    // driver.fill_region(&ferro1_region, mat::MBPgram::new(-0.3899, 0.3900, 310_000.0));
+    // driver.fill_region(&ferro2_region, mat::MBPgram::new(-0.3899, 0.3900, 310_000.0));
     // sr-latch-9; dead spot from B=[-0.03, 0.03]. This will help us see if the math is H-triggered
     // or B-triggered
-    // driver.fill_region(&ferro1_region, mat::MBFerromagnet::new(-0.3300, 0.3900, 310_000.0));
-    // driver.fill_region(&ferro2_region, mat::MBFerromagnet::new(-0.3300, 0.3900, 310_000.0));
+    // driver.fill_region(&ferro1_region, mat::MBPgram::new(-0.3300, 0.3900, 310_000.0));
+    // driver.fill_region(&ferro2_region, mat::MBPgram::new(-0.3300, 0.3900, 310_000.0));
     // mu_r=881.33, starting at H=25 to H=75.
     driver.fill_region(&ferro1_region, mat::MHPgram::new(25.0, 881.33, 44000.0));
     driver.fill_region(&ferro2_region, mat::MHPgram::new(25.0, 881.33, 44000.0));

crates/applications/archive/wrapped_torus.rs

@@ -45,7 +45,7 @@ fn main() {
     let sense1_region = Torus::new_xz(Meters::new(ferro1_center + ferro_major, half_height, half_depth), wire_major, wire_minor);
 
     //driver.fill_region(&ferro1_region, mat::db::linear_iron());
-    driver.fill_region(&ferro1_region, mat::MBFerromagnet::new(-0.3899, 0.3900, 310_000.0));
+    driver.fill_region(&ferro1_region, mat::MBPgram::new(-0.3899, 0.3900, 310_000.0));
     driver.fill_region(&drive1_region, mat::IsomorphicConductor::new(drive_conductivity));
     driver.fill_region(&sense1_region, mat::IsomorphicConductor::new(sense_conductivity));
 
@@ -54,7 +54,7 @@ fn main() {
     let drive2_region = Torus::new_xz(Meters::new(ferro2_center - ferro_major, half_height, half_depth), wire_major, wire_minor);
     let sense2_region = Torus::new_xz(Meters::new(ferro2_center + ferro_major, half_height, half_depth), wire_major, wire_minor);
 
-    driver.fill_region(&ferro2_region, mat::MBFerromagnet::new(-0.3899, 0.3900, 310_000.0));
+    driver.fill_region(&ferro2_region, mat::MBPgram::new(-0.3899, 0.3900, 310_000.0));
     driver.fill_region(&drive2_region, mat::IsomorphicConductor::new(drive_conductivity));
     driver.fill_region(&sense2_region, mat::IsomorphicConductor::new(sense_conductivity));
 

crates/applications/buffer_proto5/Cargo.toml

@@ -5,6 +5,7 @@ authors = ["Colin <colin@uninsane.org>"]
 edition = "2021"
 
 [dependencies]
+bincode = "1.3"  # MIT
 coremem = { path = "../../coremem" }
 log = "0.4"
 serde = "1.0"

crates/applications/buffer_proto5/src/cache.rs

@@ -9,6 +9,7 @@ pub struct DiskCache<K, V, S=NoSupplier> {
 }
 
 impl<K: DeserializeOwned, V: DeserializeOwned> DiskCache<K, V, NoSupplier> {
+    #[allow(dead_code)]
     pub fn new(path: &str) -> Self {
         Self::new_with_supplier(path, NoSupplier)
     }
@@ -48,6 +49,7 @@ impl<K: Serialize, V: Serialize, S> DiskCache<K, V, S> {
 }
 
 impl<K: PartialEq + Serialize, V: Serialize + Clone, S> DiskCache<K, V, S> {
+    #[allow(dead_code)]
     pub fn get_or_insert_with<F: FnOnce() -> V>(&mut self, k: K, f: F) -> V {
         if let Some(v) = self.get(&k) {
             return v.clone();

crates/applications/buffer_proto5/src/main.rs

@@ -2,23 +2,27 @@
 //! to couple them. i parameterize the entire setup over a bunch of different factors in order to
 //! search for the conditions which maximize energy transfer from the one core to the other.
 
-use coremem::{Driver, mat, meas, SpirvDriver};
-use coremem::geom::{region, Cube, Dilate, Memoize, Meters, Region, Spiral, SwapYZ, Torus, Translate, Wrap};
+use coremem::{Driver, mat, meas};
+use coremem::geom::{region, Cube, Dilate, Memoize, Meters, Spiral, SwapYZ, Torus, Translate, Wrap};
 use coremem::mat::{Ferroxcube3R1MH, IsoConductorOr};
 use coremem::real::{R32, Real as _};
 use coremem::render::CsvRenderer;
+use coremem::sim::spirv::{SpirvSim, WgpuBackend};
+use coremem::sim::units::{Seconds, Time as _};
 use coremem::stim::{CurlStimulus, Exp1, Gated, Sinusoid1, TimeVarying as _};
-use coremem::sim::units::{Seconds, Frame, Time as _};
-use coremem::sim::spirv;
-use coremem::util::cache::DiskCache;
 use log::{error, info, warn};
 use serde::{Deserialize, Serialize};
 
+mod cache;
+
+use cache::DiskCache;
+
 type Mat = IsoConductorOr<f32, Ferroxcube3R1MH>;
 
 #[allow(unused)]
 use coremem::geom::{Coord as _, Region as _};
 
+#[allow(unused)]
 #[derive(Debug, Copy, Clone, PartialEq, Eq)]
 enum PulseType {
     Square,
@@ -27,6 +31,7 @@ enum PulseType {
     ExpDecay2x,
 }
 
+#[allow(unused)]
 /// Return just the extrema of some collection
 fn extrema(mut meas: Vec<f32>) -> Vec<f32> {
     let mut i = 0;
@@ -151,24 +156,43 @@ struct Geometries {
     wrap2_len: f32,
 }
 
+/// computed measurements which get written to disk for later, manual (or grep-based) analysis.
+/// because we only write these (except for the Debug impl reading them to write to disk),
+/// rustc thinks all the fields are dead.
 #[derive(Clone, Debug, Default)]
 struct Results {
+    #[allow(dead_code)]
     m1_peak: f32,
+    #[allow(dead_code)]
     m2_peak: f32,
+    #[allow(dead_code)]
     m1_stable: f32,
+    #[allow(dead_code)]
     m2_stable: f32,
+    #[allow(dead_code)]
     h1_peak: f32,
+    #[allow(dead_code)]
     h2_max: f32,
+    #[allow(dead_code)]
     h2_min: f32,
+    #[allow(dead_code)]
     h1_stable: f32,
+    #[allow(dead_code)]
     h2_stable: f32,
+    #[allow(dead_code)]
     iset_min: f32,
+    #[allow(dead_code)]
     iset_max: f32,
+    #[allow(dead_code)]
     icoupling_peak: f32,
+    #[allow(dead_code)]
     peak_m_ratio: f32,
+    #[allow(dead_code)]
     stable_m_ratio: f32,
     /// m2_stable divided by m1_peak. i.e. "amplification"
+    #[allow(dead_code)]
     m2_stable_m1_peak: f32,
+    #[allow(dead_code)]
     t: f32,
 }
 
@@ -371,7 +395,9 @@ fn run_sim(id: u32, p: Params, g: Geometries) -> Results {
         p.clock_type,
     );
 
-    let mut driver: SpirvDriver<Mat> = Driver::new_spirv(g.dim, p.geom.feat_size);
+    let mut driver = Driver::new(SpirvSim::<f32, Mat, WgpuBackend>::new(
+        g.dim.to_index(p.geom.feat_size), p.geom.feat_size
+    ));
     driver.set_steps_per_stim(1000);
     if !driver.add_state_file(&*format!("{}/state.bc", prefix), 16000) {
         // mu_r=881.33, starting at H=25 to H=75.
@@ -396,43 +422,38 @@ fn run_sim(id: u32, p: Params, g: Geometries) -> Results {
         info!("loaded state file: skipping geometry calculations");
     }
 
-    let add_drive_sine_pulse = |driver: &mut SpirvDriver<Mat>, region: &Torus, start: f32, duration: f32, amp: f32| {
+    let add_drive_sine_pulse = |driver: &mut Driver<_>, region: &Torus, start: f32, duration: f32, amp: f32| {
         let wave = Sinusoid1::from_wavelength(amp, duration * 2.0)
             .half_cycle()
             .shifted(start);
         driver.add_stimulus(CurlStimulus::new(
             region.clone(),
             wave.clone(),
-            region.center(),
-            region.axis()
         ));
     };
 
-    let add_drive_square_pulse = |driver: &mut SpirvDriver<Mat>, region: &Torus, start: f32, duration: f32, amp: f32| {
+    let add_drive_square_pulse = |driver: &mut Driver<_>, region: &Torus, start: f32, duration: f32, amp: f32| {
         let wave = Gated::new(amp, start, start+duration);
         driver.add_stimulus(CurlStimulus::new(
             region.clone(),
             wave.clone(),
-            region.center(),
-            region.axis()
         ));
     };
 
-    let add_drive_exp_pulse = |driver: &mut SpirvDriver<Mat>, region: &Torus, start: f32, duration: f32, amp: f32| {
+    let add_drive_exp_pulse = |driver: &mut Driver<_>, region: &Torus, start: f32, duration: f32, amp: f32| {
         let wave = Exp1::new_at(amp, start, 0.5*duration);
         driver.add_stimulus(CurlStimulus::new(
             region.clone(),
             wave.clone(),
-            region.center(),
-            region.axis()
         ));
     };
 
-    let add_drive_step = |driver: &mut SpirvDriver<Mat>, region: &Torus, start: f32, amp: f32| {
-        add_drive_square_pulse(driver, region, start, 1.0, amp);
+    // step function: "permanently" increase the current by `amp`.
+    let _add_drive_step = |driver: &mut Driver<_>, region: &Torus, start: f32, amp: f32| {
+        add_drive_square_pulse(driver, region, start, 1.0 /* effectively infinite duration */, amp);
     };
 
-    let add_drive_pulse = |ty: PulseType, driver: &mut SpirvDriver<Mat>, region: &Torus, start: f32, duration: f32, amp: f32| {
+    let add_drive_pulse = |ty: PulseType, driver: &mut Driver<_>, region: &Torus, start: f32, duration: f32, amp: f32| {
         match ty {
             PulseType::Square => add_drive_square_pulse(driver, region, start, duration, amp),
             PulseType::Sine => add_drive_sine_pulse(driver, region, start, duration, amp),

crates/applications/multi_core_inverter/src/main.rs

@@ -36,21 +36,64 @@
 //! $ pushd crates/coremem; cargo run --release --bin viewer ../../out/applications/multi_core_inverter/0/ ; popd
 //! ```
 
-use coremem::geom::{Meters, Torus};
-use coremem::sim::units::Seconds;
+use coremem::geom::{Coord as _, Meters, Torus};
 use coremem::mat::{Ferroxcube3R1MH, IsoConductorOr, IsomorphicConductor};
-use coremem::spirv;
-use coremem::{Driver, SpirvDriver};
+use coremem::meas;
+use coremem::real::{self, Real as _};
+use coremem::sim::spirv::{self, SpirvSim};
+use coremem::sim::units::Seconds;
+use coremem::stim::{CurlStimulus, Exp, Gated, Sinusoid1, TimeVarying as _};
+use coremem::Driver;
 
-type Mat = IsoConductorOr<f32, Ferroxcube3R1MH>;
+type R = real::R32;
+type Mat = IsoConductorOr<R, Ferroxcube3R1MH>;
+// type Backend = spirv::CpuBackend;
+type Backend = spirv::WgpuBackend;
+type Sim = SpirvSim::<R, Mat, Backend>;
+
+#[derive(Clone, Copy)]
+enum DriveType {
+    HoldHigh,
+    ReleaseHigh,
+    HoldLow,
+    ReleaseLow,
+    Float,
+}
+
+impl DriveType {
+    fn direction(&self) -> i32 {
+        use DriveType::*;
+        match self {
+            HoldHigh | ReleaseHigh => 1,
+            HoldLow | ReleaseLow => -1,
+            Float => 0,
+        }
+    }
+    fn is_hold(&self) -> bool {
+        match self {
+            Self::HoldHigh | Self::HoldLow => true,
+            _ => false,
+        }
+    }
+    fn is_release(&self) -> bool {
+        match self {
+            Self::ReleaseHigh | Self::ReleaseLow => true,
+            _ => false,
+        }
+    }
+}
 
 fn main() {
     coremem::init_logging();
-    coremem::init_debug();
+    // coremem::init_debug();
     let um = |n| n as f32 * 1e-6;
-    let ns = |n| n as f32 * 1e-9;
+    // let ns = |n| n as f32 * 1e-9;
+    let ps = |n| n as f32 * 1e-12;
     let feat_size = um(10);
 
+    let input_magnitude = 1.0e9;
+    let clock_phase_duration = ps(2000);
+
     let s_major = um(160);
     let s_minor = um(30);
     let io_major = um(80);
@@ -61,40 +104,168 @@ fn main() {
     let sx = |n| um((n+1) * 400);
     let sy = um(400);
     let sz = um(280);
-    let couplingx = |n| sx(n) - um(200);
+    let couplingx = |n| sx(n) + um(200);
 
     let sim_bounds = Meters::new(sx(4), sy * 2.0, sz * 2.0);
     let sim_padding = Meters::new(um(80), um(80), um(80));
-    let duration = Seconds(ns(1));
 
     let drive0 = Torus::new_xz(Meters::new(sx(0) - s_major, sy, sz), io_major, io_minor);
     let sense3 = Torus::new_xz(Meters::new(sx(3) + s_major, sy, sz), io_major, io_minor);
     let ctl = |n| Torus::new_yz(Meters::new(sx(n), sy + s_major, sz), io_major, io_minor);
     let s = |n| Torus::new_xy(Meters::new(sx(n), sy, sz), s_major, s_minor);
+    // coupling(n) is the wire which couples core n into core n+1
     let coupling = |n| Torus::new_xz(Meters::new(couplingx(n), sy, sz), coupling_major, coupling_minor);
 
-    let wire_mat = IsomorphicConductor::new(1e6);
-    let ferro_mat = wire_mat;
-    // let ferro_mat = Ferroxcube3R1MH::new(); // uncomment when ready to simulate for real
+    let control_signal = |driver: &mut Driver<Sim>, region: &Torus, cycle: u32, ty: DriveType| {
+        let area = region.cross_section();
+        // negation because blog confuses current directions.
+        let amp = -ty.direction() as f32 * input_magnitude / area;
+        let start = clock_phase_duration * cycle as f32;
+        if ty.is_hold() {
+            // simple square wave:
+            let wave = Gated::new(amp, start, start + clock_phase_duration);
+            let stim = CurlStimulus::new(
+                region.clone(),
+                wave.clone(),
+            );
+            driver.add_stimulus(stim);
+        } else if ty.is_release() {
+            // decaying exponential wave:
+            let wave = Exp::new_at(amp, start, 0.125 * clock_phase_duration /* half-life */);
+            let stim = CurlStimulus::new(
+                region.clone(),
+                wave.clone(),
+            );
+            driver.add_stimulus(stim);
+        } // else: Float
+        // sine wave (untested):
+        // let wave = Sinusoid1::from_wavelength(direction as f32 * input_magnitude / area, clock_phase_duration * 2.0)
+        //     .half_cycle()
+        //     .shifted(clock_phase_duration * cycle as f32);
+    };
 
-    let mut driver: SpirvDriver<Mat> = Driver::new_spirv(sim_bounds, feat_size);
-    driver.add_classical_boundary(sim_padding);
+    let wire_mat = IsomorphicConductor::new(1e6f32.cast::<R>());
+    // let ferro_mat = wire_mat;
+    let ferro_mat = Ferroxcube3R1MH::new();
+
+    let mut driver = Driver::new(Sim::new(
+        sim_bounds.to_index(feat_size), feat_size,
+    ));
+    driver.add_classical_boundary_explicit::<R, _>(sim_padding);
+
+    //////// create the wires and toroids
     driver.fill_region(&drive0, wire_mat);
     driver.fill_region(&sense3, wire_mat);
     for core in 0..4 {
         driver.fill_region(&s(core), ferro_mat);
         driver.fill_region(&ctl(core), wire_mat);
-        if core != 0 {
+        if core != 3 {
             driver.fill_region(&coupling(core), wire_mat);
         }
     }
 
-    let prefix = "out/applications/multi_core_inverter/0/";
+    //////// monitor some measurements
+    // driver.add_measurement(meas::CurrentLoop::new("drv0", drive0.clone()));
+    for core in 0..4 {
+        driver.add_measurement(meas::CurrentLoop::new(
+            &format!("drive{}", core),
+            ctl(core),
+        ));
+    }
+    for core in 0..3 {
+        driver.add_measurement(meas::CurrentLoop::new(
+            &format!("sense{}", core),
+            coupling(core),
+        ));
+    }
+    for core in 0..4 {
+        driver.add_measurement(meas::MagneticLoop::new(
+            &format!("state{}", core),
+            s(core),
+        ));
+    }
+    driver.add_measurement(meas::CurrentLoop::new("sense3", sense3.clone()));
+
+    //////// create the stimuli
+
+    // each row N denotes the drive currents at clock cycle N.
+    // each col M denotes the drive current at core M.
+    // this drive map is taken from the blog article as of 2022-08-09; it has noted errors
+    // let drive_map = [
+    //     [1, 0, 1, 1],
+    //     [1, 1, 0, 1],
+    //     [1, 1, 1, 0],
+    //     [0, 1, 1, 1i32],
+    // ];
+    // proposed drive map: the location where a core is charged to 1 requires a negative current.
+    // let drive_map = [
+    //     [ 1,  0,  1,  1],    // clear S0 -> S1; hold S2, S3
+    //     [-1,  1,  0,  1],    // clear S1 -> S2; hold S3; reset S0
+    //     [-1,  1,  1,  0],    // clear S2 -> S3; hold S0 (high), S1
+    //     [ 0,  1, -1,  1i32], // clear S3 -> Sn; hold S1, S2 (S0 gets new input
+    // ];
+    // the above drive map has some implicit requirement of clock lag. it won't work with perfectly
+    // square clocks. here's a realizable drive map:
+    // #[allow(non_snake_case)]
+    // let F = 0; // indicates a gate which has been driven to a known value and is now floating
+    // #[allow(non_snake_case)]
+    // let D = 0; // indicates a gate which carries data.
+    // let drive_map = [
+    //     [-1,  1,  1,  1],    // initial state: S0=1, S1=S2=S3=0
+    //     [-1,  F,  1,  1],    // open S1 for write
+
+    //     [ 1,  D,  1,  1],    // clear S0 -> S1; hold S2, S3
+    //     [ 1,  D,  F,  1],    // open S2 for write
+
+    //     [-1,  1,  D,  1],    // clear S1 -> S2; hold S3; reset S0
+    //     [-1,  1,  D,  F],    // open S3 for write
+
+    //     [-1,  1,  1,  D],    // clear S2 -> S3; hold S0 (high), S1
+    //     [ F,  1,  1,  D],    // open S0 for write (not used in this demo)
+
+    //     [ D,  1,  1,  1i32], // clear S3 -> Sn; hold S1, S2 (S0 gets new input)
+    //     [ D,  F,  1,  1i32], // open S1 for write (not used; for completeness)
+    // ];
+
+    use DriveType::*;
+    let drive_map = [
+        [HoldLow,     HoldHigh,     HoldHigh,     HoldHigh],    // initial state: S0=1, S1=S2=S3=0
+        [HoldLow,     ReleaseHigh,  HoldHigh,     HoldHigh],    // open S1 for write
+
+        [HoldHigh,    Float,        HoldHigh,     HoldHigh],    // clear S0 -> S1; hold S2, S3
+        [HoldHigh,    Float,        ReleaseHigh,  HoldHigh],    // open S2 for write
+
+        [HoldHigh,    HoldHigh,     Float,        HoldHigh],    // clear S1 -> S2; hold S0, S3
+        // NB: bonus clock! we have to reset S0 before opening S3, else writeback from next stage
+        [HoldLow,     HoldHigh,     Float,        HoldHigh],    // reset S0; hold S1, S3
+        [HoldLow,     HoldHigh,     Float,        ReleaseHigh], // open S3 for write
+
+        [HoldLow,     HoldHigh,     HoldHigh,     Float],       // clear S2 -> S3; hold S0 (high), S1
+        [ReleaseLow,  HoldHigh,     HoldHigh,     Float],       // open S0 for write (not used in this demo)
+
+        [Float,       HoldHigh,     HoldHigh,     HoldHigh],    // clear S3 -> Sn; hold S1, S2 (S0 gets new input)
+        [Float,       ReleaseHigh,  HoldHigh,     HoldHigh],    // open S1 for write (not used; for completeness)
+    ];
+
+    let cycles = drive_map.len() as u32;
+    let duration = Seconds(clock_phase_duration * (cycles + 3) as f32);
+
+    for cycle in 0..cycles {
+        for core in 0..4 {
+            // current polarity of the drive wires is inverted in this model v.s. the blog article,
+            // so invert our currents in order to achieve the magnetic states of the article.
+            let sig = drive_map[cycle as usize][core as usize];
+            control_signal(&mut driver, &ctl(core), cycle, sig);
+        }
+    }
+
+    let prefix = "out/applications/multi_core_inverter/8/";
     let _ = std::fs::create_dir_all(&prefix);
     // driver.add_state_file(&*format!("{}state.bc", prefix), 9600);
-    driver.add_serializer_renderer(&*format!("{}frame-", prefix), 36000, None);
-    driver.add_csv_renderer(&*format!("{}meas.csv", prefix), 400, None);
-    driver.set_steps_per_stim(200);
+    driver.add_serializer_renderer(&*format!("{}frame-", prefix), 400, None);
+    driver.add_csv_renderer(&*format!("{}meas-detailed.csv", prefix), 100, None);
+    driver.add_csv_renderer(&*format!("{}meas.csv", prefix), 800, None);
+    driver.set_steps_per_stim(100);
 
     driver.step_until(duration);
 }

crates/applications/sr_latch/src/main.rs

@@ -3,9 +3,9 @@
 /// the SR latch in this example is wired to a downstream latch, mostly to show that it's
 /// possible to transfer the state (with some limitation) from one latch to another.
 
-use coremem::{Driver, mat, meas, SpirvDriver};
-use coremem::geom::{Meters, Torus};
-use coremem::sim::spirv;
+use coremem::{Driver, mat, meas};
+use coremem::geom::{Coord as _, Meters, Torus};
+use coremem::sim::spirv::{SpirvSim, WgpuBackend};
 use coremem::sim::units::Seconds;
 use coremem::stim::{CurlStimulus, Sinusoid1, TimeVarying as _};
 
@@ -59,7 +59,9 @@ fn main() {
     let coupling_region = Torus::new_xz(Meters::new(0.5*(ferro1_center + ferro2_center), ferro_center_y, half_depth), wire_coupling_major, wire_minor);
     let sense_region = Torus::new_xz(Meters::new(ferro2_center + ferro_major, ferro_center_y, half_depth), wire_major, wire_minor);
 
-    let mut driver: SpirvDriver<mat::FullyGenericMaterial<f32>> = Driver::new_spirv(Meters::new(width, height, depth), feat_size);
+    let mut driver = Driver::new(SpirvSim::<f32, mat::FullyGenericMaterial<f32>, WgpuBackend>::new(
+        Meters::new(width, height, depth).to_index(feat_size), feat_size
+    ));
 
     // mu_r=881.33, starting at H=25 to H=75.
     driver.fill_region(&ferro1_region, mat::MHPgram::new(25.0, 881.33, 44000.0));
@@ -82,8 +84,6 @@ fn main() {
         driver.add_stimulus(CurlStimulus::new(
             region.clone(),
             wave.clone(),
-            region.center(),
-            region.axis()
         ));
     };
 

crates/applications/wavefront/src/main.rs

@@ -7,10 +7,14 @@
 //! with something that absorbs energy. since this example doesn't, it lets you see what
 //! happens when you just use the default boundary conditions.
 
-use coremem::{mat, driver};
-use coremem::geom::{Coord as _, Cube, Index, Vec3};
+use coremem::{mat, Driver};
+use coremem::geom::{Coord as _, Cube, Index};
 use coremem::units::Seconds;
+use coremem::sim::spirv::{self, SpirvSim};
 use coremem::stim::{Stimulus, TimeVarying as _, UniformStimulus};
+use coremem::cross::vec::Vec3;
+
+type Mat = mat::FullyGenericMaterial<f32>;
 
 fn main() {
     coremem::init_logging();
@@ -21,14 +25,13 @@ fn main() {
     // each cell represents 1um x 1um x 1um volume
     let feature_size = 1e-6;
 
-    // Create the simulation "driver" which uses the CPU as backend.
-    // by default all the computations are done with R32: a f32 which panics on NaN/Inf
-    // you can parameterize it to use R64, or unchecked f32 -- see src/driver.rs for the definition
-    let mut driver: driver::CpuDriver = driver::Driver::new(size, feature_size);
-
-    // uncomment to use the Spirv/GPU driver. this one is restricted to unchecked f32.
-    // note: this won't have better perf unless you reduce the y4m/term renderer framerate below.
-    // let mut driver: driver::SpirvDriver = driver::Driver::new_spirv(size, feature_size);
+    // create the simulation "driver".
+    // the first parameter is the float type to use: f32 for unchecked math, coremem::real::R32
+    // to guard against NaN/Inf (useful for debugging).
+    // to run this on the gpu instead of the gpu, replace `CpuBackend` with `WgpuBackend`.
+    let mut driver = Driver::new(SpirvSim::<f32, Mat, spirv::CpuBackend>::new(
+        size, feature_size
+    ));
 
     // create a conductor on the left side.
     let conductor = Cube::new(

crates/coremem/Cargo.toml

@@ -16,7 +16,6 @@ crossterm = "0.24"  # MIT
 csv = "1.1"  # MIT or Unlicense
 dashmap = "5.3"  # MIT
 dyn-clone = "1.0"  # MIT or Apache 2.0
-enum_dispatch = "0.3"  # MIT or Apache 2.0
 env_logger = "0.9"  # MIT or Apache 2.0
 float_eq = "1.0"  # MIT or Apache 2.0
 font8x8 = "0.3"  # MIT
@@ -47,7 +46,7 @@ spirv-std = { git = "https://github.com/EmbarkStudios/rust-gpu" }
 spirv-std-macros = { git = "https://github.com/EmbarkStudios/rust-gpu" }
 spirv_backend = { path = "../spirv_backend" }
 spirv_backend_runner = { path = "../spirv_backend_runner" }
-coremem_types = { path = "../types", features = ["fmt", "serde"] }
+coremem_cross = { path = "../cross", features = ["fmt", "serde"] }
 
 
 [dev-dependencies]

crates/coremem/benches/driver.rs

@@ -1,27 +1,16 @@
-use coremem::{Driver, SimState, SpirvDriver};
+use coremem::Driver;
 use coremem::geom::Index;
-use coremem::mat::{Ferroxcube3R1MH, IsoConductorOr, GenericMaterial, GenericMaterialNoPml, GenericMaterialOneField};
-use coremem::real::R32;
-use coremem::sim::spirv::{self, SpirvSim};
+use coremem::mat::{Ferroxcube3R1MH, IsoConductorOr, FullyGenericMaterial};
+use coremem::sim::spirv::{SpirvSim, WgpuBackend};
 use criterion::{BenchmarkId, criterion_group, criterion_main, Criterion};
 
-type DefaultDriver = Driver::<SimState<R32, GenericMaterial<R32>>>;
-
-pub fn bench_step(c: &mut Criterion) {
-    for size in &[10, 20, 40, 80, 160] {
-        let sim = SimState::<R32, GenericMaterial<R32>>::new(Index::new(*size, *size, *size), 1e-5);
-        c.bench_with_input(BenchmarkId::new("Driver::step", size), &sim, |b, sim| {
-            let mut driver = Driver::new_with_state(sim.clone());
-            b.iter(|| driver.step())
-        });
-    }
-}
 
 pub fn bench_step_spirv(c: &mut Criterion) {
+    type Mat = FullyGenericMaterial<f32>;
     for size in &[10, 20, 40, 80, 160] {
-        let sim: SpirvSim = SpirvSim::new(Index::new(*size, *size, *size), 1e-5);
+        let sim = SpirvSim::<f32, Mat, WgpuBackend>::new(Index::new(*size, *size, *size), 1e-5);
         c.bench_with_input(BenchmarkId::new("Driver::step_spirv", size), &sim, |b, sim| {
-            let mut driver = Driver::new_with_state(sim.clone());
+            let mut driver = Driver::new(sim.clone());
             b.iter(|| driver.step())
         });
     }
@@ -30,42 +19,13 @@ pub fn bench_step_spirv(c: &mut Criterion) {
 pub fn bench_step_spirv_iso_3r1(c: &mut Criterion) {
     type Mat = IsoConductorOr<f32, Ferroxcube3R1MH>;
     for size in &[10, 20, 40, 80, 160] {
-        let sim: SpirvSim<Mat> = SpirvSim::new(Index::new(*size, *size, *size), 1e-5);
+        let sim = SpirvSim::<f32, Mat, WgpuBackend>::new(Index::new(*size, *size, *size), 1e-5);
         c.bench_with_input(BenchmarkId::new("Driver::spirv_ISO3R1", size), &sim, |b, sim| {
-            let mut driver: SpirvDriver<Mat> = Driver::new_with_state(sim.clone());
+            let mut driver = Driver::new(sim.clone());
             b.iter(|| driver.step())
         });
     }
 }
 
-// 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 = DefaultDriver::new(Index::new(size, size, size), 1e-5);
-//             b.iter(|| driver.step())
-//         });
-//     }
-// }
-// 
-// 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 = DefaultDriver::new(Index::new(size, size, size), 1e-5);
-//             b.iter(|| driver.step())
-//         });
-//     }
-// }
-// 
-// 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 = DefaultDriver::new(Index::new(size, size, size), 1e-5);
-//             driver.add_pml_boundary(Index::new(thickness, thickness, thickness));
-//             b.iter(|| driver.step())
-//         });
-//     }
-// }
-
-criterion_group!(benches, /*bench_step,*/ bench_step_spirv, bench_step_spirv_iso_3r1);
+criterion_group!(benches, bench_step_spirv, bench_step_spirv_iso_3r1);
 criterion_main!(benches);

crates/coremem/benches/history.txt

@@ -103,3 +103,15 @@ Driver::spirv_ISO3R1/20 time:   [548.70 us 555.85 us 563.28 us]
 Driver::spirv_ISO3R1/40 time:   [1.5333 ms 1.5405 ms 1.5489 ms]
 Driver::spirv_ISO3R1/80 time:   [13.299 ms 13.335 ms 13.376 ms]
 Driver::spirv_ISO3R1/160time:   [164.57 ms 164.74 ms 164.93 ms]
+
+5a0766451d96835061a674ab94f00341adb2b187:
+Driver::step_spirv/10   time:   [590.26 us 600.42 us 613.28 us]
+Driver::step_spirv/20   time:   [870.49 us 884.81 us 902.21 us]
+Driver::step_spirv/40   time:   [3.4094 ms 3.4285 ms 3.4498 ms]
+Driver::step_spirv/80   time:   [35.488 ms 35.673 ms 35.922 ms]
+Driver::step_spirv/160  time:   [270.98 ms 271.19 ms 271.43 ms]
+Driver::spirv_ISO3R1/10 time:   [585.57 us 596.11 us 608.79 us]
+Driver::spirv_ISO3R1/20 time:   [826.63 us 841.79 us 860.86 us]
+Driver::spirv_ISO3R1/40 time:   [2.8808 ms 2.9004 ms 2.9237 ms]
+Driver::spirv_ISO3R1/80 time:   [28.955 ms 29.027 ms 29.115 ms]
+Driver::spirv_ISO3R1/160time:   [216.03 ms 216.22 ms 216.45 ms]

crates/coremem/src/bin/bench.rs

@@ -1,5 +1,7 @@
-use coremem::{self, Driver, GenericSim, SimState};
-use coremem::sim::spirv::SpirvSim;
+use coremem::{self, Driver, AbstractSim};
+use coremem::sim::legacy::SimState;
+use coremem::sim::spirv::{SpirvSim, WgpuBackend};
+use coremem::cross::mat::FullyGenericMaterial;
 use coremem::geom::Index;
 use std::time::{Instant, Duration};
 
@@ -16,19 +18,35 @@ fn measure<F: FnMut()>(name: &str, n_times: u32, mut f: F) -> f32 {
     avg
 }
 
-fn measure_steps<S: GenericSim + Clone + Default + Send + Sync + 'static>(name: &str, steps_per_call: u32, mut d: Driver<S>) {
+fn measure_steps<S: AbstractSim + Clone + Default + Send + 'static>(name: &str, steps_per_call: u32, mut d: Driver<S>) {
     measure(name, 100/steps_per_call, || d.step_multiple(steps_per_call));
 }
 
 fn main() {
     coremem::init_logging();
-    measure_steps("spirv/80", 1, Driver::<SpirvSim>::new_spirv(Index::new(80, 80, 80), 1e-3));
-    measure_steps("sim/80", 1, Driver::<SimState>::new(Index::new(80, 80, 80), 1e-3));
-    measure_steps("spirv/80 step(2)", 2, Driver::<SpirvSim>::new_spirv(Index::new(80, 80, 80), 1e-3));
-    measure_steps("sim/80 step(2)", 2, Driver::<SimState>::new(Index::new(80, 80, 80), 1e-3));
-    measure_steps("spirv/80 step(10)", 10, Driver::<SpirvSim>::new_spirv(Index::new(80, 80, 80), 1e-3));
-    measure_steps("sim/80 step(10)", 10, Driver::<SimState>::new(Index::new(80, 80, 80), 1e-3));
-    measure_steps("spirv/80 step(100)", 100, Driver::<SpirvSim>::new_spirv(Index::new(80, 80, 80), 1e-3));
-    measure_steps("sim/80 step(100)", 100, Driver::<SimState>::new(Index::new(80, 80, 80), 1e-3));
+    measure_steps("spirv/80", 1, Driver::new(
+        SpirvSim::<f32, FullyGenericMaterial<f32>, WgpuBackend>::new(Index::new(80, 80, 80), 1e-3)
+    ));
+    measure_steps("sim/80", 1, Driver::<SimState>::new(
+        SimState::new(Index::new(80, 80, 80), 1e-3)
+    ));
+    measure_steps("spirv/80 step(2)", 2, Driver::new(
+        SpirvSim::<f32, FullyGenericMaterial<f32>, WgpuBackend>::new(Index::new(80, 80, 80), 1e-3)
+    ));
+    measure_steps("sim/80 step(2)", 2, Driver::<SimState>::new(
+        SimState::new(Index::new(80, 80, 80), 1e-3)
+    ));
+    measure_steps("spirv/80 step(10)", 10, Driver::new(
+        SpirvSim::<f32, FullyGenericMaterial<f32>, WgpuBackend>::new(Index::new(80, 80, 80), 1e-3)
+    ));
+    measure_steps("sim/80 step(10)", 10, Driver::<SimState>::new(
+        SimState::new(Index::new(80, 80, 80), 1e-3)
+    ));
+    measure_steps("spirv/80 step(100)", 100, Driver::new(
+        SpirvSim::<f32, FullyGenericMaterial<f32>, WgpuBackend>::new(Index::new(80, 80, 80), 1e-3)
+    ));
+    measure_steps("sim/80 step(100)", 100, Driver::<SimState>::new(
+        SimState::new(Index::new(80, 80, 80), 1e-3)
+    ));
 }
 

crates/coremem/src/bin/pml_tuning.rs

@@ -1,496 +0,0 @@
-use coremem::*;
-use coremem::geom::*;
-use coremem::real::{Real as _, ToFloat as _};
-use coremem::stim::AbstractStimulus;
-use rand::rngs::StdRng;
-use rand::{Rng as _, SeedableRng as _};
-
-fn energy<R: Region>(s: &dyn SampleableSim, reg: &R) -> f32 {
-    let e = f64::half() * s.map_sum_over_enumerated(reg, |_pos: Index, cell| {
-        cell.e().mag_sq().to_f64()
-    });
-    e.cast()
-}
-
-fn energy_now_and_then<R: Region>(state: &mut StaticSim, reg: &R, frames: u32) -> (f32, f32) {
-    let energy_0 = energy(state, reg);
-    for _ in 0..frames {
-        state.step();
-    }
-    let energy_1 = energy(state, reg);
-    (energy_0, energy_1)
-}
-
-struct PmlStim<F> {
-    /// Maps index -> (stim vector, stim frequency)
-    f: F,
-    t_end: f32,
-    feat_size: f32,
-}
-impl<F: Fn(Index) -> (Vec3<f32>, f32) + Sync> AbstractStimulus for PmlStim<F> {
-    fn at(&self, t_sec: f32, pos: Meters) -> (Vec3<f32>, Vec3<f32>) {
-        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;
-        let sig = sig_angle.sin();
-        (e*gate*sig, Vec3::zero())
-    }
-}
-
-/// Apply some stimulus, and then let it decay and measure the ratio of energy left in the system
-fn apply_stim_full_interior<F>(state: &mut StaticSim, frames: u32, f: F)
-where F: Fn(Index) -> (Vec3<f32>, f32) + Sync // returns (E vector, omega)
-{
-    let stim = PmlStim {
-        f,
-        t_end: (frames as f32) * state.timestep(),
-        feat_size: state.feature_size(),
-    };
-    for _t in 0..frames {
-        state.apply_stimulus(&stim);
-        state.step();
-    }
-}
-
-fn apply_stim_over_region<R, F>(state: &mut StaticSim, frames: u32, reg: R, f: F)
-where
-    R: Region,
-    F: Fn(Index) -> (Vec3<f32>, f32) + Sync,
-{
-    let feat = state.feature_size();
-    apply_stim_full_interior(state, frames, |idx| {
-        if reg.contains(idx.to_meters(feat)) {
-            f(idx)
-        } else {
-            (Vec3::zero(), 0.0)
-        }
-    });
-}
-
-/// Stimulate each point in the region with a pseudorandom (but predictable) e wave
-fn apply_chaotic_stim_over_region<R: Region>(state: &mut StaticSim, frames: u32, interior: R) {
-    apply_stim_over_region(state, frames, interior, |idx| {
-        let seed = (idx.x() as u64) ^ ((idx.y() as u64) << 16) ^ ((idx.z() as u64) << 32);
-        let mut rng = StdRng::seed_from_u64(seed);
-        let dir = Vec3::new(
-            rng.gen_range(-1.0..1.0),
-            rng.gen_range(-1.0..1.0),
-            rng.gen_range(-1.0..1.0),
-        );
-        // XXX only works if it's a whole number. I suppose this makes sense though, as
-        // other numbers would create higher harmonics when gated.
-        let hz = rng.gen_range(0..=2);
-        (dir, hz as _)
-    })
-}
-
-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;
-        let lower_right_idx = state.size() - Index::unit() - upper_left_idx;
-        let interior = Cube::new(upper_left_idx.to_meters(feat), lower_right_idx.to_meters(feat));
-        apply_chaotic_stim_over_region(state, frames, interior);
-    }
-
-    let upper_left_idx = Index::unit()*boundary;
-    let lower_right_idx = state.size() - Index::unit() - upper_left_idx;
-    let sim_region = Cube::new(upper_left_idx.to_meters(feat), lower_right_idx.to_meters(feat));
-
-    let (energy_0, energy_1) = energy_now_and_then(state, &sim_region, frames);
-    // println!("Energy: {}/{}", energy_1, energy_0);
-    energy_1/energy_0
-}
-
-#[allow(unused)]
-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(pow);
-        let coord_stretch = b * sc_coeff / timestep;
-        let conductivity = Vec3::unit() * (b.mag() * cond_coeff / timestep);
-        unimplemented!();
-        Static::default()
-        // Static {
-        //     // TODO PML coord_stretch,
-        //     conductivity,
-        //     pml: Some((PmlState::new(), PmlParameters::new(coord_stretch))),
-        //     ..Default::default()
-        // }
-    });
-    state
-}
-
-fn main() {
-    // Explanation here (slide 63): https://empossible.net/wp-content/uploads/2020/01/Lecture-The-Perfectly-Matched-Layer.pdf
-    // Claims that eps0/delta_t * n^3 is a good way to activate the PML layer
-
-    for pow in vec![3.0] {
-    for sc_coeff in vec![
-        // 0.1*consts::EPS0,
-        // 0.5*consts::EPS0,
-        // 1e0*consts::EPS0,
-        // 1e1*consts::EPS0,
-        // 1e2*consts::EPS0,
-        // 1e3*consts::EPS0,
-        // 1e6*consts::EPS0,
-        // 1e9*consts::EPS0,
-        // 1e12*consts::EPS0,
-        // 1e15*consts::EPS0,
-        // 1e18*consts::EPS0,
-        // 1e21*consts::EPS0,
-        // 0.01,
-        // 0.03,
-        // 0.05,
-        // 0.07,
-        // 0.08,
-        // 0.09,
-        // 0.1,
-        // 0.15,
-        // 0.2,
-        // 0.25,
-        // 0.3,
-        // 0.4,
-        // 0.5,
-        // 0.6,
-        // 0.7,
-        // 0.8,
-        // 0.9,
-        //0.95,
-        // 1.0,
-        // 1.5,
-        // 2.0,
-        // 3.0,
-        // 5.0,
-        // 10.0,
-        // 100.0,
-        // 1000.0,
-        //0.0,
-        0.5,
-    ] {
-    //for cond_coeff in vec![0.0, 1.0, 1e3, 1e6, 1e9] {
-    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] {
-                    let size = Index::unit()*121;
-                    let sim_inset = 40;
-                    let boundary = Index::unit()*pml_width;
-                    let mut state = state_for_pml(size, boundary, feat_size, sc_coeff, cond_coeff, pow);
-                    let ratio = chaotic_pml_test(&mut state, pml_width, sim_inset, frames);
-                    println!("{},{} (pow={}, f={}, width={}, frames={}): {}", sc_coeff, cond_coeff, pow, feat_size, pml_width, frames, ratio);
-                }
-            }
-        }
-    }
-    }
-    }
-    // Conclusions:
-    // * if coordinate stretching is divided by time_step, then the absorption is independent of
-    //   feature size.
-    // For 240 frames, PML width of 20, size=121, sim_inset = 40, cubic onset:
-    // * coefficients between [0.4, 1.0] are all within 30% of each other
-    // * coefficients between [0.5, 1.0] are all within 20% of each other
-    // * coefficients > 1 show instability (i.e. instable if coord_stretch > 1 / timestep);
-    // * absorption generally increases as the coefficient approaches 1.0 from the left.
-    //   This begins to reverse at least by 0.98
-    // For 2400 frames, instability starts somewhere between 0.7 and 0.8
-    //   0.5 (f=0.000001, width=20, frames=2400): 0.0006807010986857421
-    //   0.6 (f=0.000001, width=20, frames=2400): 0.0005800974138738364
-    //   0.7 (f=0.000001, width=20, frames=2400): 0.0005092274390086559
-    //   0.8 (f=0.000001, width=20, frames=2400): 7326609898580109000000
-    //   These numbers use hz = rng.gen_range(0..=5), with 20 frames of excitation and 20 steps
-    //     between source and boundary
-    //
-    // Reducing hz to 0..=2 achieves much better perf (pow=3.0):
-    //   0.5 (f=0.000001, width=20, frames=2400): 0.0000007135657380376141
-    //   0.6 (f=0.000001, width=20, frames=2400): 0.0000006561934691721446
-    //   0.7 (f=0.000001, width=20, frames=2400): 0.0000006099334150762627
-    //   0.8 (f=0.000001, width=20, frames=2400): 2180868728529433400000
-    // Quadratic (pow=2.0):
-    //   0.4 (pow=2, f=0.000001, width=20, frames=2400): 0.0000006839243836328471
-    //   0.5 (pow=2, f=0.000001, width=20, frames=2400): 0.0000006189051555210391
-    //   0.6 (pow=2, f=0.000001, width=20, frames=2400): 0.9331062414894028
-    //   0.7 (pow=2, f=0.000001, width=20, frames=2400): 222818414359827930000000000000000000000000000000000000000000000000000000000000000000000
-    // Possibly better efficacy, but also less stable
-    // Linear (pow=1.0) is 100% unstable for at least coeff >= 0.5
-    // Subcubic (pow=2.5):
-    //   0.4 (pow=2.5, f=0.000001, width=20, frames=2400): 0.0000007333781495449756
-    //   0.5 (pow=2.5, f=0.000001, width=20, frames=2400): 0.0000006632381361190291
-    //   0.6 (pow=2.5, f=0.000001, width=20, frames=2400): 0.0000006089982278171636
-    //   0.7 (pow=2.5, f=0.000001, width=20, frames=2400): 7853952559.189004
-    // Superlinear (pow=1.5):
-    //   0.4 (pow=1.5, f=0.000001, width=20, frames=2400): 0.0000006502259131444893
-    //   0.5 (pow=1.5, f=0.000001, width=20, frames=2400): 0.0000005956833662247809
-    //   0.6 (pow=1.5, f=0.000001, width=20, frames=2400): 175615904583581400000000000000000000000000000000000000000000000000000000000000000000
-    // pow=4.0:
-    //   0.4 (pow=4, f=0.000001, width=20, frames=2400): 0.0000009105207758718423
-    //   0.5 (pow=4, f=0.000001, width=20, frames=2400): 0.0000008169331995553951
-    //   0.6 (pow=4, f=0.000001, width=20, frames=2400): 0.0000007525813705681166
-    //   0.7 (pow=4, f=0.000001, width=20, frames=2400): 0.0000007019565296008103
-    //   0.8 (pow=4, f=0.000001, width=20, frames=2400): 0.0000006600546205898854
-    //   0.9 (pow=4, f=0.000001, width=20, frames=2400): 0.0000006246047577819345
-    //   1 (pow=4, f=0.000001, width=20, frames=2400): 5437370254206590000000000000000000000000000000000
-    // pow=3.5:
-    //   0.4 (pow=3.5, f=0.000001, width=20, frames=2400): 0.0000008485707012112478
-    //   0.5 (pow=3.5, f=0.000001, width=20, frames=2400): 0.0000007653501316913873
-    //   0.6 (pow=3.5, f=0.000001, width=20, frames=2400): 0.0000007047188963430957
-    //   0.7 (pow=3.5, f=0.000001, width=20, frames=2400): 0.0000006561795324449849
-    //   0.8 (pow=3.5, f=0.000001, width=20, frames=2400): 0.0000006159621438553542
-    //   0.9 (pow=3.5, f=0.000001, width=20, frames=2400): 18913062838424785000000000000000000
-    // So generally lower powers = more absorbtion (makes sense: higher average coordinate
-    // stretching), but is more sensitive to error. All powers > 1.5 are stable at coeff=0.5.
-    // We don't know that much about reflection from just this data.
-    // Running over a smaller timeframe should give some suggestion about reflection
-    //   0.4 (pow=2, f=0.000001, width=20, frames=120): 0.2870021971493659
-    //   0.5 (pow=2, f=0.000001, width=20, frames=120): 0.2647114916517679 **
-    //   0.6 (pow=2, f=0.000001, width=20, frames=120): 0.24807341265362437
-    //   0.7 (pow=2, f=0.000001, width=20, frames=120): 0.2350081813369852
-    //   0.8 (pow=2, f=0.000001, width=20, frames=120): 0.22438565227661014
-    //   0.9 (pow=2, f=0.000001, width=20, frames=120): 0.21552516261496907
-    //   1.0 (pow=2, f=0.000001, width=20, frames=120): 0.20798702383197107
-    //   1.5 (pow=2, f=0.000001, width=20, frames=120): 0.39444643124947426
-    //   0.4 (pow=3, f=0.000001, width=20, frames=120): 0.3596522778473902
-    //   0.5 (pow=3, f=0.000001, width=20, frames=120): 0.33718409856439474
-    //   0.6 (pow=3, f=0.000001, width=20, frames=120): 0.32025414395335816
-    //   0.7 (pow=3, f=0.000001, width=20, frames=120): 0.3067653188201421 **
-    //   0.8 (pow=3, f=0.000001, width=20, frames=120): 0.2956250635204507
-    //   0.9 (pow=3, f=0.000001, width=20, frames=120): 0.2861895432992242
-    //   1.0 (pow=3, f=0.000001, width=20, frames=120): 0.27804573152917056
-    //   1.5 (pow=3, f=0.000001, width=20, frames=120): 0.24950413937515897
-    //   0.4 (pow=4, f=0.000001, width=20, frames=120): 0.41043160705678194
-    //   0.5 (pow=4, f=0.000001, width=20, frames=120): 0.38844802581052806
-    //   0.6 (pow=4, f=0.000001, width=20, frames=120): 0.3721037027266112
-    //   0.7 (pow=4, f=0.000001, width=20, frames=120): 0.35912262947037576
-    //   0.8 (pow=4, f=0.000001, width=20, frames=120): 0.34837556188585944
-    //   0.9 (pow=4, f=0.000001, width=20, frames=120): 0.33922665697709947 **
-    //   1 (pow=4, f=0.000001, width=20, frames=120): 0.33128127577349487
-    //   1.5 (pow=4, f=0.000001, width=20, frames=120): 0.30260541581006584
-    // Even smaller timeframe:
-    //   0.5 (pow=2, f=0.000001, width=20, frames=5): 1.0109467232603757
-    //   0.5 (pow=2, f=0.000001, width=20, frames=10): 1.0117049984675424
-    //   0.5 (pow=2, f=0.000001, width=20, frames=20): 1.015054295037722
-    //   0.5 (pow=2, f=0.000001, width=20, frames=30): 1.016127315217072
-    //   0.5 (pow=2, f=0.000001, width=20, frames=40): 1.0125450428623481
-    //   0.5 (pow=2, f=0.000001, width=20, frames=60): 0.9520114341681021
-    //   0.5 (pow=2, f=0.000001, width=20, frames=80): 0.7617289087518861
-    //   0.5 (pow=2, f=0.000001, width=20, frames=120): 0.2647114916517679
-    //   0.5 (pow=2, f=0.000001, width=20, frames=160): 0.018698432021826004
-    //   0.5 (pow=2, f=0.000001, width=20, frames=200): 0.001053437888276037
-    //   0.5 (pow=3, f=0.000001, width=20, frames=5): 1.0109467232603757
-    //   0.5 (pow=3, f=0.000001, width=20, frames=10): 1.0117049984674478
-    //   0.5 (pow=3, f=0.000001, width=20, frames=20): 1.0150542142468075
-    //   0.5 (pow=3, f=0.000001, width=20, frames=30): 1.0161755433699602
-    //   0.5 (pow=3, f=0.000001, width=20, frames=40): 1.0140118995944478
-    //   0.5 (pow=3, f=0.000001, width=20, frames=60): 0.9828724159361404
-    //   0.5 (pow=3, f=0.000001, width=20, frames=80): 0.8304056107875255
-    //   0.5 (pow=3, f=0.000001, width=20, frames=120): 0.33718409856439474
-    //   0.5 (pow=3, f=0.000001, width=20, frames=160): 0.033906877503334515
-    //   0.5 (pow=3, f=0.000001, width=20, frames=200): 0.0016075121270576818
-    //   0.5 (pow=4, f=0.000001, width=20, frames=5): 1.0109467232603757
-    //   0.5 (pow=4, f=0.000001, width=20, frames=10): 1.0117049984674424
-    //   0.5 (pow=4, f=0.000001, width=20, frames=20): 1.015054203569572
-    //   0.5 (pow=4, f=0.000001, width=20, frames=30): 1.0161827617273314
-    //   0.5 (pow=4, f=0.000001, width=20, frames=40): 1.014355067000562
-    //   0.5 (pow=4, f=0.000001, width=20, frames=60): 0.9971421502814346
-    //   0.5 (pow=4, f=0.000001, width=20, frames=80): 0.8718457443603458
-    //   0.5 (pow=4, f=0.000001, width=20, frames=120): 0.38844802581052806
-    //   0.5 (pow=4, f=0.000001, width=20, frames=160): 0.04879003869310861
-    //   0.5 (pow=4, f=0.000001, width=20, frames=200): 0.0025039595751290534
-    //   That the numbers are all the same for t < 20 is not encouraging.
-    //     Could be because of the courant number 0.577? No energy has reached the border yet?
-    // Why don't we query just the energy in the sim region -- not the boundary?
-    // After filtering to measure energy only in the sim region (note that some energy outside the
-    // sim region COULD be reflected back into the sim region later):
-    //   0.5 (pow=2, f=0.000001, width=20, frames=5): 1.0109467232603757
-    //   0.5 (pow=2, f=0.000001, width=20, frames=10): 1.0117049984565956
-    //   0.5 (pow=2, f=0.000001, width=20, frames=20): 1.0150051760172862
-    //   0.5 (pow=2, f=0.000001, width=20, frames=30): 1.0103422463150569
-    //   0.5 (pow=2, f=0.000001, width=20, frames=40): 0.9534926741875023
-    //   0.5 (pow=2, f=0.000001, width=20, frames=60): 0.7052236436305864
-    //   0.5 (pow=2, f=0.000001, width=20, frames=80): 0.42366103692252166
-    //   0.5 (pow=2, f=0.000001, width=20, frames=120): 0.047117482623349645
-    //   0.5 (pow=2, f=0.000001, width=20, frames=160): 0.0012788501490854916
-    //   0.5 (pow=2, f=0.000001, width=20, frames=200): 0.00025052555946608536
-    //   0.5 (pow=2, f=0.000001, width=20, frames=300): 0.000026163786700826282
-    //   0.5 (pow=2, f=0.000001, width=20, frames=400): 0.000009482696335385068
-    //   0.5 (pow=2, f=0.000001, width=20, frames=600): 0.0000034339203626681873
-    //   0.5 (pow=2, f=0.000001, width=20, frames=800): 0.0000016501593549063537
-    //   0.5 (pow=2, f=0.000001, width=20, frames=1200): 0.0000007097133927819365
-    //   0.5 (pow=2, f=0.000001, width=20, frames=1600): 0.0000004990355277480898
-    //   0.5 (pow=2, f=0.000001, width=20, frames=2400): 0.00000026541631012862064
-    //   0.5 (pow=3, f=0.000001, width=20, frames=5): 1.0109467232603757
-    //   0.5 (pow=3, f=0.000001, width=20, frames=10): 1.0117049984564943
-    //   0.5 (pow=3, f=0.000001, width=20, frames=20): 1.015005103609574
-    //   0.5 (pow=3, f=0.000001, width=20, frames=30): 1.010347173332254
-    //   0.5 (pow=3, f=0.000001, width=20, frames=40): 0.9534945122673749
-    //   0.5 (pow=3, f=0.000001, width=20, frames=60): 0.7052084303225945
-    //   0.5 (pow=3, f=0.000001, width=20, frames=80): 0.42365590263359737
-    //   0.5 (pow=3, f=0.000001, width=20, frames=120): 0.04711930602404218
-    //   0.5 (pow=3, f=0.000001, width=20, frames=160): 0.0013058690825149086
-    //   0.5 (pow=3, f=0.000001, width=20, frames=200): 0.00027124952640186995
-    //   0.5 (pow=3, f=0.000001, width=20, frames=300): 0.00003363606115185493
-    //   0.5 (pow=3, f=0.000001, width=20, frames=400): 0.000012979216745112146
-    //   0.5 (pow=3, f=0.000001, width=20, frames=600): 0.000004790890253991059
-    //   0.5 (pow=3, f=0.000001, width=20, frames=800): 0.0000018967555372300478
-    //   0.5 (pow=3, f=0.000001, width=20, frames=1200): 0.000000785458387440694
-    //   0.5 (pow=3, f=0.000001, width=20, frames=1600): 0.0000005657588895614973
-    //   0.5 (pow=3, f=0.000001, width=20, frames=2400): 0.0000002931443936827707
-    //   0.5 (pow=4, f=0.000001, width=20, frames=5): 1.0109467232603757
-    //   0.5 (pow=4, f=0.000001, width=20, frames=10): 1.011704998456489
-    //   0.5 (pow=4, f=0.000001, width=20, frames=20): 1.0150050988039427
-    //   0.5 (pow=4, f=0.000001, width=20, frames=30): 1.0103477261980267
-    //   0.5 (pow=4, f=0.000001, width=20, frames=40): 0.9534933477754826
-    //   0.5 (pow=4, f=0.000001, width=20, frames=60): 0.7052093991486525
-    //   0.5 (pow=4, f=0.000001, width=20, frames=80): 0.4236543609610083
-    //   0.5 (pow=4, f=0.000001, width=20, frames=120): 0.047143352151935734
-    //   0.5 (pow=4, f=0.000001, width=20, frames=160): 0.0014253338264054373
-    //   0.5 (pow=4, f=0.000001, width=20, frames=200): 0.0003719236701949888
-    //   0.5 (pow=4, f=0.000001, width=20, frames=300): 0.0000655709065994179
-    //   0.5 (pow=4, f=0.000001, width=20, frames=400): 0.000019428436023699366
-    //   0.5 (pow=4, f=0.000001, width=20, frames=600): 0.000006787101433521501
-    //   0.5 (pow=4, f=0.000001, width=20, frames=800): 0.000002352253311002263
-    //   0.5 (pow=4, f=0.000001, width=20, frames=1200): 0.0000008945009475225972
-    //   0.5 (pow=4, f=0.000001, width=20, frames=1600): 0.0000006298015155592009
-    //   0.5 (pow=4, f=0.000001, width=20, frames=2400): 0.0000003110877898243663
-    //   pow=2 and pow=3 look nearly interchangeable in perf, with pow=4 slightly worse.
-    //   Given higher stability for pow=3, it pushes me in that direction
-    // What about ordinary conductivity?
-    //   Uniaxial conductors do nothing
-    // non-axial conductors (measured over full volume -- not just sim volume)
-    //   0,0.000000000000001 (pow=2, f=0.000001, width=20, frames=20): 1.0150542017357447
-    //   0,0.000000000000001 (pow=2, f=0.000001, width=20, frames=40): 1.0144968836020802
-    //   0,0.000000000000001 (pow=2, f=0.000001, width=20, frames=80): 1.001249693707569
-    //   0,0.000000000000001 (pow=2, f=0.000001, width=20, frames=120): 0.9703627391410313
-    //   0,0.000000000000001 (pow=2, f=0.000001, width=20, frames=240): 0.9839347439326915
-    //   0,0.000000000000001 (pow=2, f=0.000001, width=20, frames=400): 0.9710971333858839
-    //   0,0.000000000001 (pow=2, f=0.000001, width=20, frames=20): 1.015054152058065
-    //   0,0.000000000001 (pow=2, f=0.000001, width=20, frames=40): 1.0139965206685178
-    //   0,0.000000000001 (pow=2, f=0.000001, width=20, frames=80): 0.8611234324272511
-    //   0,0.000000000001 (pow=2, f=0.000001, width=20, frames=120): 0.38724876605299674
-    //   0,0.000000000001 (pow=2, f=0.000001, width=20, frames=240): 0.010452374686360136
-    //   0,0.000000000001 (pow=2, f=0.000001, width=20, frames=400): 0.0006325928742141608
-    //   0,0.000000001 (pow=2, f=0.000001, width=20, frames=20): 1.0150333781336882
-    //   0,0.000000001 (pow=2, f=0.000001, width=20, frames=40): 0.9834360538048873
-    //   0,0.000000001 (pow=2, f=0.000001, width=20, frames=80): 0.5072413415421974
-    //   0,0.000000001 (pow=2, f=0.000001, width=20, frames=120): 0.10345998021237998
-    //   0,0.000000001 (pow=2, f=0.000001, width=20, frames=240): 0.007519161562708659
-    //   0,0.000000001 (pow=2, f=0.000001, width=20, frames=400): 0.0006070229440826888
-    //   0,0.000001 (pow=2, f=0.000001, width=20, frames=20): 1.0150056763216164
-    //   0,0.000001 (pow=2, f=0.000001, width=20, frames=40): 1.0123907067781637
-    //   0,0.000001 (pow=2, f=0.000001, width=20, frames=80): 1.0095377475794036
-    //   0,0.000001 (pow=2, f=0.000001, width=20, frames=120): 1.0027226600640793
-    //   0,0.000001 (pow=2, f=0.000001, width=20, frames=240): 1.002880617201309
-    //   0,0.000001 (pow=2, f=0.000001, width=20, frames=400): 0.9817776859704778
-    //   0,0.001 (pow=2, f=0.000001, width=20, frames=20): 1.0150056260987463
-    //   0,0.001 (pow=2, f=0.000001, width=20, frames=40): 1.0127202910529511
-    //   0,0.001 (pow=2, f=0.000001, width=20, frames=80): 1.0141889645709163
-    //   0,0.001 (pow=2, f=0.000001, width=20, frames=120): 1.012712629869813
-    //   0,0.001 (pow=2, f=0.000001, width=20, frames=240): 1.0225762750160987
-    //   0,0.001 (pow=2, f=0.000001, width=20, frames=400): 1.0158169138074233
-    //   0,1 (pow=2, f=0.000001, width=20, frames=20): 1.0150056260488272
-    //   0,1 (pow=2, f=0.000001, width=20, frames=40): 1.0127206223077903
-    //   0,1 (pow=2, f=0.000001, width=20, frames=80): 1.01419363999915
-    //   0,1 (pow=2, f=0.000001, width=20, frames=120): 1.0127226953239867
-    //   0,1 (pow=2, f=0.000001, width=20, frames=240): 1.0225963209694455
-    //   0,1 (pow=2, f=0.000001, width=20, frames=400): 1.0158520285873147
-    //   0,0.000000000000001 (pow=3, f=0.000001, width=20, frames=20): 1.0150542017810857
-    //   0,0.000000000000001 (pow=3, f=0.000001, width=20, frames=40): 1.014497284230248
-    //   0,0.000000000000001 (pow=3, f=0.000001, width=20, frames=80): 1.0013349322815948
-    //   0,0.000000000000001 (pow=3, f=0.000001, width=20, frames=120): 0.9707653838346726
-    //   0,0.000000000000001 (pow=3, f=0.000001, width=20, frames=240): 0.98542720202683
-    //   0,0.000000000000001 (pow=3, f=0.000001, width=20, frames=400): 0.9736668760490382
-    //   0,0.000000000001 (pow=3, f=0.000001, width=20, frames=20): 1.0150541972999592
-    //   0,0.000000000001 (pow=3, f=0.000001, width=20, frames=40): 1.0143834319294165
-    //   0,0.000000000001 (pow=3, f=0.000001, width=20, frames=80): 0.9087702844329181
-    //   0,0.000000000001 (pow=3, f=0.000001, width=20, frames=120): 0.4641021861485614
-    //   0,0.000000000001 (pow=3, f=0.000001, width=20, frames=240): 0.025542565246607578
-    //   0,0.000000000001 (pow=3, f=0.000001, width=20, frames=400): 0.0024678463432975324
-    //   0,0.000000001 (pow=3, f=0.000001, width=20, frames=20): 1.0150507306076173
-    //   0,0.000000001 (pow=3, f=0.000001, width=20, frames=40): 0.9973368290451762
-    //   0,0.000000001 (pow=3, f=0.000001, width=20, frames=80): 0.5359561255382069
-    //   0,0.000000001 (pow=3, f=0.000001, width=20, frames=120): 0.10107229161449789
-    //   0,0.000000001 (pow=3, f=0.000001, width=20, frames=240): 0.0018726822600725027
-    //   0,0.000000001 (pow=3, f=0.000001, width=20, frames=400): 0.0001395247365920271
-    //   0,0.000001 (pow=3, f=0.000001, width=20, frames=20): 1.015006719517991
-    //   0,0.000001 (pow=3, f=0.000001, width=20, frames=40): 1.0067208901320701
-    //   0,0.000001 (pow=3, f=0.000001, width=20, frames=80): 0.9295573543881088
-    //   0,0.000001 (pow=3, f=0.000001, width=20, frames=120): 0.8372562924822498
-    //   0,0.000001 (pow=3, f=0.000001, width=20, frames=240): 0.7109091146305759
-    //   0,0.000001 (pow=3, f=0.000001, width=20, frames=400): 0.5405933200528408
-    //   0,0.001 (pow=3, f=0.000001, width=20, frames=20): 1.015005627048278
-    //   0,0.001 (pow=3, f=0.000001, width=20, frames=40): 1.0127139915480314
-    //   0,0.001 (pow=3, f=0.000001, width=20, frames=80): 1.0141000517606538
-    //   0,0.001 (pow=3, f=0.000001, width=20, frames=120): 1.0125212238134818
-    //   0,0.001 (pow=3, f=0.000001, width=20, frames=240): 1.0221951620616994
-    //   0,0.001 (pow=3, f=0.000001, width=20, frames=400): 1.0151495300965783
-    //   0,1 (pow=3, f=0.000001, width=20, frames=20): 1.0150056260497768
-    //   0,1 (pow=3, f=0.000001, width=20, frames=40): 1.012720616007617
-    //   0,1 (pow=3, f=0.000001, width=20, frames=80): 1.0141935510766382
-    //   0,1 (pow=3, f=0.000001, width=20, frames=120): 1.0127225038875007
-    //   0,1 (pow=3, f=0.000001, width=20, frames=240): 1.0225959397081388
-    //   0,1 (pow=3, f=0.000001, width=20, frames=400): 1.0158513607157937
-    //   It appears here that the cubic roll-off is best, and very low conductivities are required.
-    //
-    // Isotropic conductor, energy measured only within the sim area:
-    //   The optimized case (0.5*EPS0/timestep * x^3) is almost IDENTICAL to the
-    //   optimized stretched-coordinate version.
-    //   0,0.000000000000001 (pow=2, f=0.000001, width=20, frames=120): 0.07433102240945513
-    //   0,0.000000000000001 (pow=2, f=0.000001, width=20, frames=400): 0.20431918328907037
-    //   0,0.000000000001 (pow=2, f=0.000001, width=20, frames=120): 0.04851013262210713
-    //   0,0.000000000001 (pow=2, f=0.000001, width=20, frames=400): 0.00018366232272528987
-    //   0,0.0000000000044270939064065 (pow=2, f=0.000001, width=20, frames=120): 0.047101055930619994
-    //   0,0.0000000000044270939064065 (pow=2, f=0.000001, width=20, frames=400): 0.000007399479477711689
-    //   0,0.000000000008854187812813 (pow=2, f=0.000001, width=20, frames=120): 0.047132676278627536
-    //   0,0.000000000008854187812813 (pow=2, f=0.000001, width=20, frames=400): 0.000008105989741898135
-    //   0,0.00000000001 (pow=2, f=0.000001, width=20, frames=120): 0.047145523858393434
-    //   0,0.00000000001 (pow=2, f=0.000001, width=20, frames=400): 0.000008327266269591764
-    //   0,0.0000000001 (pow=2, f=0.000001, width=20, frames=120): 0.05005311352870215
-    //   0,0.0000000001 (pow=2, f=0.000001, width=20, frames=400): 0.00003665711687010239
-    //   0,0.000000001 (pow=2, f=0.000001, width=20, frames=120): 0.0808038989604745
-    //   0,0.000000001 (pow=2, f=0.000001, width=20, frames=400): 0.0005636469832522345
-    //   0,0.00000001 (pow=2, f=0.000001, width=20, frames=120): 0.42090062694288544
-    //   0,0.00000001 (pow=2, f=0.000001, width=20, frames=400): 0.07558396270665715
-    //   0,0.000001 (pow=2, f=0.000001, width=20, frames=120): 0.9456301791158588
-    //   0,0.000001 (pow=2, f=0.000001, width=20, frames=400): 0.9456217562930543
-    //   0,0.001 (pow=2, f=0.000001, width=20, frames=120): 0.9549600788853474
-    //   0,0.001 (pow=2, f=0.000001, width=20, frames=400): 0.9782784298162882
-    //   0,1 (pow=2, f=0.000001, width=20, frames=120): 0.9549694776242036
-    //   0,1 (pow=2, f=0.000001, width=20, frames=400): 0.9783121245194134
-    //   0,0.000000000000001 (pow=3, f=0.000001, width=20, frames=120): 0.07435108171392021
-    //   0,0.000000000000001 (pow=3, f=0.000001, width=20, frames=400): 0.20479144240774633
-    //   0,0.000000000001 (pow=3, f=0.000001, width=20, frames=120): 0.050005777221700895
-    //   0,0.000000000001 (pow=3, f=0.000001, width=20, frames=400): 0.0006548443359188922
-    //   0,0.0000000000044270939064065 (pow=3, f=0.000001, width=20, frames=120): 0.0471097820852083
-    //   0,0.0000000000044270939064065 (pow=3, f=0.000001, width=20, frames=400): 0.00000726683693079002
-    //   0,0.000000000008854187812813 (pow=3, f=0.000001, width=20, frames=120): 0.04711795151894901
-    //   0,0.000000000008854187812813 (pow=3, f=0.000001, width=20, frames=400): 0.000007527865420649209
-    //   0,0.00000000001 (pow=3, f=0.000001, width=20, frames=120): 0.0471214615521427
-    //   0,0.00000000001 (pow=3, f=0.000001, width=20, frames=400): 0.000007561968840324101
-    //   0,0.0000000001 (pow=3, f=0.000001, width=20, frames=120): 0.04775494018209187
-    //   0,0.0000000001 (pow=3, f=0.000001, width=20, frames=400): 0.00001655968038405813
-    //   0,0.000000001 (pow=3, f=0.000001, width=20, frames=120): 0.05548225210925576
-    //   0,0.000000001 (pow=3, f=0.000001, width=20, frames=400): 0.00011973516700797696
-    //   0,0.00000001 (pow=3, f=0.000001, width=20, frames=120): 0.09859866928724702
-    //   0,0.00000001 (pow=3, f=0.000001, width=20, frames=400): 0.0010788775694009376
-    //   0,0.000001 (pow=3, f=0.000001, width=20, frames=120): 0.7905518369602526
-    //   0,0.000001 (pow=3, f=0.000001, width=20, frames=400): 0.5217962235754597
-    //   0,0.001 (pow=3, f=0.000001, width=20, frames=120): 0.9547813505519078
-    //   0,0.001 (pow=3, f=0.000001, width=20, frames=400): 0.9776380394998598
-    //   0,1 (pow=3, f=0.000001, width=20, frames=120): 0.9549692988681137
-    //   0,1 (pow=3, f=0.000001, width=20, frames=400): 0.978311483657099
-    // With both PML and conductor boundary:
-    //   0.5,0.0000000000044270939064065 (pow=2, f=0.000001, width=20, frames=120): 0.04709454626708049
-    //   0.5,0.0000000000044270939064065 (pow=2, f=0.000001, width=20, frames=400): 0.000007626399231684735
-    //   0.5,0.0000000000044270939064065 (pow=2, f=0.000001, width=20, frames=1200): 0.0000008116632120088743
-    //   0.5,0.0000000000044270939064065 (pow=2, f=0.000001, width=20, frames=2400): 0.00000032497678753300923
-    //   0.5,0.0000000000044270939064065 (pow=3, f=0.000001, width=20, frames=120): 0.04709779630205149
-    //   0.5,0.0000000000044270939064065 (pow=3, f=0.000001, width=20, frames=400): 0.000007132448872463483
-    //   0.5,0.0000000000044270939064065 (pow=3, f=0.000001, width=20, frames=1200): 0.000000710451109532798
-    //   0.5,0.0000000000044270939064065 (pow=3, f=0.000001, width=20, frames=2400): 0.0000002515070714149805
-    //   This is basically no change from JUST PML or JUST conductors
-    //   Maybe I should be trying to vary the width: maybe PML works more effectively for narrower
-    //   boundaries than do conductors?
-}

crates/coremem/src/driver.rs

@@ -1,12 +1,11 @@
-use crate::geom::{Coord, Index, Meters, Region, Vec3};
-use crate::mat::{self, Pml};
+use crate::geom::{Coord, Index, Meters, Region};
+use crate::mat;
 use crate::meas::{self, AbstractMeasurement};
-use crate::real::{self, Real};
+use crate::real::Real;
 use crate::render::{self, MultiRenderer, Renderer};
-use crate::sim::{GenericSim, MaterialSim, SampleableSim, SimState};
+use crate::sim::AbstractSim;
 use crate::sim::units::{Frame, Time};
-use crate::sim::spirv::{self, SpirvSim};
-use crate::stim::AbstractStimulus;
+use crate::stim::{self, AbstractStimulus};
 
 use log::{info, trace};
 use serde::{Deserialize, Serialize};
@@ -16,18 +15,19 @@ use std::sync::mpsc::{sync_channel, SyncSender, Receiver};
 use std::time::{Duration, Instant};
 use threadpool::ThreadPool;
 
-pub struct Driver<S=SimState> {
-    pub state: S,
+pub struct Driver<S> {
+    state: S,
     renderer: Arc<MultiRenderer<S>>,
     // TODO: use Rayon's thread pool?
     render_pool: ThreadPool,
     render_channel: (SyncSender<()>, Receiver<()>),
+    // TODO: split out into a Diagnostics struct
     time_spent_stepping: Duration,
     time_spent_on_stimuli: Duration,
     time_spent_prepping_render: Duration,
     time_spent_blocked_on_render: Duration,
     time_spent_rendering: Arc<Mutex<Duration>>,
-    measurements: Vec<Box<dyn AbstractMeasurement>>,
+    measurements: Vec<Arc<dyn AbstractMeasurement<S>>>,
     stimuli: StimuliAdapter,
     start_time: Instant,
     last_diag_time: Instant,
@@ -35,25 +35,8 @@ pub struct Driver<S=SimState> {
     sim_end_time: Option<Frame>,
 }
 
-pub type CpuDriver<R=real::R32, M=mat::GenericMaterial<R>> = Driver<SimState<R, M>>;
-pub type SpirvDriver<M=mat::FullyGenericMaterial<f32>> = Driver<SpirvSim<M>>;
-
-impl<R: Real, M: Default> Driver<SimState<R, M>> {
-    pub fn new<C: Coord>(size: C, feature_size: f32) -> Self {
-        Self::new_with_state(SimState::new(size.to_index(feature_size), feature_size))
-    }
-}
-
-impl<M: spirv::IntoFfi> SpirvDriver<M>
-    where M::Ffi: Default + 'static
-{
-    pub fn new_spirv<C: Coord>(size: C, feature_size: f32) -> Self {
-        Self::new_with_state(SpirvSim::new(size.to_index(feature_size), feature_size))
-    }
-}
-
-impl<S> Driver<S> {
-    pub fn new_with_state(state: S) -> Self {
+impl<S: AbstractSim> Driver<S> {
+    pub fn new(state: S) -> Self {
         Self {
             state,
             renderer: Arc::new(MultiRenderer::new()),
@@ -65,10 +48,10 @@ impl<S> Driver<S> {
             time_spent_blocked_on_render: Default::default(),
             time_spent_rendering: Default::default(),
             measurements: vec![
-                Box::new(meas::Time),
-                Box::new(meas::Meta),
-                Box::new(meas::Energy::world()),
-                Box::new(meas::Power::world()),
+                Arc::new(meas::Time),
+                Arc::new(meas::Meta),
+                Arc::new(meas::Energy::world()),
+                Arc::new(meas::Power::world()),
             ],
             stimuli: StimuliAdapter::new(),
             start_time: Instant::now(),
@@ -76,15 +59,13 @@ impl<S> Driver<S> {
             sim_end_time: None,
         }
     }
-}
 
-impl<S> Driver<S> {
     pub fn add_stimulus<Stim: AbstractStimulus + 'static>(&mut self, s: Stim) {
         self.stimuli.push(Box::new(s))
     }
 
-    pub fn add_measurement<Meas: AbstractMeasurement + 'static>(&mut self, m: Meas) {
-        self.measurements.push(Box::new(m));
+    pub fn add_measurement<Meas: AbstractMeasurement<S> + 'static>(&mut self, m: Meas) {
+        self.measurements.push(Arc::new(m));
     }
 
     pub fn set_steps_per_stim(&mut self, steps_per_stim: u64) {
@@ -92,16 +73,20 @@ impl<S> Driver<S> {
     }
 }
 
-impl<S: MaterialSim> Driver<S> {
+impl<S: AbstractSim> Driver<S> {
     pub fn fill_region<Reg: Region, M: Into<S::Material> + Clone>(&mut self, region: &Reg, mat: M) {
         self.state.fill_region(region, mat);
     }
-    pub fn test_region_filled<Reg: Region, M: Into<S::Material> + Clone>(&mut self, region: &Reg, mat: M) -> bool {
+    pub fn test_region_filled<Reg: Region, M>(&mut self, region: &Reg, mat: M) -> bool
+    where
+        M: Into<S::Material> + Clone,
+        S::Material: PartialEq
+    {
         self.state.test_region_filled(region, mat)
     }
 }
 
-impl<S: SampleableSim> Driver<S> {
+impl<S: AbstractSim> Driver<S> {
     pub fn size(&self) -> Index {
         self.state.size()
     }
@@ -113,11 +98,7 @@ impl<S: SampleableSim> Driver<S> {
     }
 }
 
-impl<S: SampleableSim + Send + Sync + 'static> Driver<S> {
-    pub fn dyn_state(&mut self) -> &mut dyn SampleableSim {
-        &mut self.state
-    }
-
+impl<S: AbstractSim + 'static> Driver<S> {
     fn add_renderer<Rend: Renderer<S> + 'static>(
         &mut self, renderer: Rend, name: &str, step_frequency: u64, frame_limit: Option<u64>
     ) {
@@ -140,31 +121,27 @@ impl<S: SampleableSim + Send + Sync + 'static> Driver<S> {
     }
 }
 
-impl<S: SampleableSim + Send + Sync + Serialize + 'static> Driver<S> {
+impl<S: AbstractSim + Serialize + 'static> Driver<S> {
     pub fn add_serializer_renderer(&mut self, out_base: &str, step_frequency: u64, frame_limit: Option<u64>) {
         let fmt_str = format!("{out_base}{{step_no}}.bc", out_base=out_base);
-        self.add_renderer(render::SerializerRenderer::new_static(&*fmt_str), &*fmt_str, step_frequency, frame_limit);
+        self.add_renderer(render::SerializerRenderer::new_generic(&*fmt_str), &*fmt_str, step_frequency, frame_limit);
     }
 }
 
-impl<S: SampleableSim + Send + Sync + Serialize + for<'a> Deserialize<'a> + 'static> Driver<S> {
+impl<S: AbstractSim + Send + Sync + Serialize + for<'a> Deserialize<'a> + 'static> Driver<S> {
     /// instruct the driver to periodically save the simulation state to the provided path.
     /// also attempts to load an existing state file, returning `true` on success.
     pub fn add_state_file(&mut self, state_file: &str, snapshot_frequency: u64) -> bool {
         let ser = render::SerializerRenderer::new(state_file);
-        let loaded = match ser.try_load() {
-            Some(state) => {
-                self.state = state.state;
-                true
-            },
-            None => false,
-        };
+        let loaded = ser.try_load().map(|s| {
+            self.state = s.state;
+        }).is_some();
         self.add_renderer(ser, state_file, snapshot_frequency, None);
         loaded
     }
 }
 
-impl<S: GenericSim + Clone + Default + Send + Sync + 'static> Driver<S> {
+impl<S: AbstractSim + Clone + Default + Send + 'static> Driver<S> {
     fn render(&mut self) {
         let prep_start = Instant::now();
         let their_state = self.state.clone();
@@ -177,7 +154,8 @@ impl<S: GenericSim + Clone + Default + Send + Sync + 'static> Driver<S> {
             sender.send(()).unwrap();
             trace!("render begin");
             let start_time = Instant::now();
-            renderer.render(&their_state, &*their_measurements, Default::default());
+            let meas: Vec<&dyn AbstractMeasurement<S>> = their_measurements.iter().map(|m| &**m).collect();
+            renderer.render(&their_state, &*meas, Default::default());
             *time_spent_rendering.lock().unwrap() += start_time.elapsed();
             trace!("render end");
         });
@@ -261,7 +239,7 @@ impl<S: GenericSim + Clone + Default + Send + Sync + 'static> Driver<S> {
         let sim_end_time = sim_end_time.to_frame(self.state.timestep());
         self.sim_end_time = Some(sim_end_time);
         let mut stepped = false;
-        while self.dyn_state().step_no() < *sim_end_time {
+        while self.state.step_no() < *sim_end_time {
             self.step_multiple(100);
             stepped = true;
         }
@@ -274,18 +252,7 @@ impl<S: GenericSim + Clone + Default + Send + Sync + 'static> Driver<S> {
     }
 }
 
-impl<S: MaterialSim> Driver<S> {
-    pub fn add_pml_boundary<C: Coord, R: Real>(&mut self, thickness: C)
-        where S::Material: From<Pml<R>>
-    {
-        let timestep = self.state.timestep();
-        self.state.fill_boundary_using(thickness, |boundary_ness| {
-            let b = boundary_ness.elem_pow(3.0);
-            let conductivity = b * (0.5 / timestep);
-            Pml::new(conductivity)
-        });
-    }
-
+impl<S: AbstractSim> Driver<S> {
     pub fn add_classical_boundary<C: Coord>(&mut self, thickness: C)
         where S::Material: From<mat::IsomorphicConductor<f32>>
     {
@@ -319,7 +286,7 @@ struct StimuliAdapter {
 }
 
 impl AbstractStimulus for StimuliAdapter {
-    fn at(&self, t_sec: f32, pos: Meters) -> (Vec3<f32>, Vec3<f32>) {
+    fn at(&self, t_sec: f32, pos: Meters) -> stim::Fields {
         self.stim.at(t_sec, pos)
         // TODO: remove this stuff  (here only for testing)
         /*

crates/coremem/src/geom/line.rs

@@ -1,5 +1,5 @@
-use crate::geom::Vec2;
-use crate::real::Real;
+use coremem_cross::real::Real;
+use coremem_cross::vec::Vec2;
 
 use std::ops::Add;
 

crates/coremem/src/geom/mod.rs

@@ -6,8 +6,22 @@ mod units;
 pub use line::Line2d;
 pub use polygon::Polygon2d;
 pub use region::{
-    Cube, CylinderZ, Dilate, InvertedRegion, Memoize, Region, Sphere, Spiral, SwapXZ, SwapYZ, Torus, Translate, Union, WorldRegion, Wrap
+    Cube,
+    CylinderZ,
+    Dilate,
+    HasCrossSection,
+    InvertedRegion,
+    Memoize,
+    Region,
+    Sphere,
+    Spiral,
+    SwapXZ,
+    SwapYZ,
+    Torus,
+    Translate,
+    Union,
+    WorldRegion,
+    Wrap,
 };
 pub use units::{Coord, Meters, OrdMeters, Index};
-pub use coremem_types::vec::{Vec2, Vec3, Vec3u};
 

crates/coremem/src/geom/polygon.rs

@@ -1,5 +1,6 @@
-use crate::geom::{Line2d, Vec2};
-use crate::real::Real;
+use crate::geom::Line2d;
+use coremem_cross::real::Real;
+use coremem_cross::vec::Vec2;
 
 #[derive(Clone, Debug, PartialEq)]
 pub struct Polygon2d<R> {

crates/coremem/src/geom/region/mod.rs

@@ -1,5 +1,7 @@
 use crate::geom::{Coord, Meters, OrdMeters};
 
+use coremem_cross::vec::Vec3;
+
 use dyn_clone::{self, DynClone};
 use rayon::prelude::*;
 use serde::{Serialize, Deserialize};
@@ -16,6 +18,13 @@ pub trait Region: Send + Sync + DynClone {
 }
 dyn_clone::clone_trait_object!(Region);
 
+/// some (volume) which has a tangent vector everywhere inside/on it.
+/// for example, a cylinder has tangents everywhere except its axis.
+/// the return vector should be normalized, or zero.
+pub trait HasCrossSection {
+    fn cross_section_normal(&self, p: Meters) -> Vec3<f32>;
+}
+
 pub fn and<T1: Region + 'static, T2: Region + 'static>(r1: T1, r2: T2) -> Intersection {
     Intersection::new().and(r1).and(r2)
 }

crates/coremem/src/geom/region/primitives.rs

@@ -1,11 +1,12 @@
-use crate::geom::{Meters, Vec2, Vec3};
-use crate::real::Real as _;
+use crate::geom::Meters;
+use coremem_cross::real::Real as _;
+use coremem_cross::vec::{Vec2, Vec3};
 
 use serde::{Serialize, Deserialize};
 use std::fmt::{self, Display};
 use std::ops::Range;
 
-use super::Region;
+use super::{HasCrossSection, Region};
 
 #[derive(Copy, Clone, Serialize, Deserialize)]
 pub struct CylinderZ {
@@ -86,7 +87,7 @@ impl Region for Torus {
         // 2. Find the point `q` on the circle which is nearest to `p`.
         // 3. Consider the distance from `p` to `q`.
         let rel_p = *p - *self.center;
-        let p_on_plane = rel_p - self.normal.with_mag(self.normal.dot(rel_p));
+        let p_on_plane = rel_p - self.normal.with_mag(self.normal.dot(rel_p)).unwrap();
         let q = if p_on_plane == Vec3::zero() {
             // avoid division by zero.
             // The point is precisely on the axis of the torus.
@@ -94,15 +95,24 @@ impl Region for Torus {
             // and they all give the same answer.
             // Such a point is given by rotating the normal axis by 90 degrees in ANY DIRECTION
             let off_axis = self.normal.arbitrary_orthogonal_vector();
-            off_axis.with_mag(self.major_rad)
+            off_axis.with_mag(self.major_rad).unwrap()
         } else {
-            p_on_plane.with_mag(self.major_rad)
+            p_on_plane.with_mag(self.major_rad).unwrap()
         };
         let distance_to_circle_sq = rel_p.distance_sq(q);
         distance_to_circle_sq < self.minor_rad * self.minor_rad
     }
 }
 
+impl HasCrossSection for Torus {
+    fn cross_section_normal(&self, coord: Meters) -> Vec3<f32> {
+        let axis = self.axis();
+        let to_coord = *coord - *self.center();
+        // this creates a normal which always points "counter-clockwise" along the shape
+        axis.cross(to_coord).norm()
+    }
+}
+
 #[derive(Copy, Clone, Serialize, Deserialize)]
 pub struct Sphere {
     center: Meters,

crates/coremem/src/geom/units.rs

@@ -1,6 +1,6 @@
-use crate::real::ToFloat;
+use coremem_cross::real::ToFloat;
+use coremem_cross::vec::{Vec3, Vec3u};
 use serde::{Serialize, Deserialize};
-use super::{Vec3, Vec3u};
 use std::fmt::{self, Display};
 use std::cmp::Ordering;
 use std::ops::{Add, Deref, Div, Mul, Neg, Sub};

crates/coremem/src/lib.rs

@@ -9,17 +9,16 @@ use log::info;
 
 pub mod driver;
 pub mod geom;
-pub mod mat;
 pub mod meas;
 pub mod render;
 pub mod sim;
 pub mod stim;
-pub mod util;
 
 pub use driver::*;
-pub use mat::*;
 pub use sim::*;
-pub use coremem_types::real;
+pub use coremem_cross as cross;
+pub use coremem_cross::real;
+pub use coremem_cross::mat;
 
 // Some things to keep in mind:
 // B = mu_r*H + M

crates/coremem/src/mat/mb_ferromagnet.rs

@@ -1,195 +0,0 @@
-use crate::geom::{Line2d, Vec2, Vec3};
-use crate::mat::Material;
-use crate::real::Real;
-use crate::sim::CellState;
-
-use serde::{Serialize, Deserialize};
-
-/// M(B) parallelogram
-///
-///```text
-///     ____________
-///    /           /
-///   /     .     /
-///  /           /
-/// /___________/
-/// ```
-///
-/// The `.` depicts (0, 0). X axis is B; y axis is M.
-/// As B increases, M remains constant until it hits an edge.
-/// Then M rises up to its max.
-/// Same thing happens on the left edge, as B decreases and M falls to its min.
-#[derive(Default, Copy, Clone, PartialEq, Serialize, Deserialize)]
-pub struct MBPgram<R> {
-    /// Vertical range of the graph
-    pub max_m: R,
-    /// X coordinate at which the upward slope starts
-    pub b_start: R,
-    /// X coordinate at which the upward slope ends
-    pub b_end: R,
-}
-
-impl<R: Real> MBPgram<R> {
-    pub fn new(b_start: R, b_end: R, max_m: R) -> Self {
-        Self { b_start, b_end, max_m }
-    }
-
-    /// Return the new `M`
-    pub fn move_b(&self, m: R, target_b: R) -> R {
-        let right_edge = Line2d::new(
-            Vec2::new(self.b_start, -self.max_m),
-            Vec2::new(self.b_end, self.max_m),
-        );
-        let left_edge = Line2d::new(
-            Vec2::new(-self.b_start, self.max_m),
-            Vec2::new(-self.b_end, -self.max_m),
-        );
-        // m must be at least this much:
-        let min_m = right_edge.clamp_by_x(target_b).y();
-        // m must be no more than this:
-        let max_m = left_edge.clamp_by_x(target_b).y();
-        m.max_or_undefined(min_m).min_or_undefined(max_m)
-    }
-}
-
-#[derive(Default, Copy, Clone, PartialEq, Serialize, Deserialize)]
-pub struct NativeMBFerromagnet<R> {
-    m: Vec3<R>,
-    curve: MBPgram<R>,
-}
-
-impl<R: Real> NativeMBFerromagnet<R> {
-    pub fn new<R2: Real>(b_start: R2, b_end: R2, max_m: R2) -> Self {
-        Self {
-            m: Vec3::zero(),
-            curve: MBPgram::new(b_start.cast(), b_end.cast(), max_m.cast()),
-        }
-    }
-
-    pub fn curve(&self) -> MBPgram<R> {
-        self.curve
-    }
-}
-
-impl<R: Real> Material<R> for NativeMBFerromagnet<R> {
-    fn step_b(&mut self, context: &CellState<R>, delta_b: Vec3<R>) {
-        let target_b = context.with_m(self.m).b() + delta_b;
-        // println!("step_b {}", target_b);
-        self.m = Vec3::new(
-            self.curve.move_b(self.m.x(), target_b.x()),
-            self.curve.move_b(self.m.y(), target_b.y()),
-            self.curve.move_b(self.m.z(), target_b.z()),
-        );
-    }
-    fn m(&self) -> Vec3<R> {
-        self.m
-    }
-}
-
-#[derive(Default, Copy, Clone, PartialEq, Serialize, Deserialize)]
-pub struct SpirvMBFerromagnet<R>(NativeMBFerromagnet<R>);
-
-impl<R: Real> SpirvMBFerromagnet<R> {
-    pub fn new<R2: Real>(b_start: R2, b_end: R2, max_m: R2) -> Self {
-        Self(NativeMBFerromagnet::new(b_start, b_end, max_m))
-    }
-
-    pub fn curve(&self) -> MBPgram<R> {
-        self.0.curve()
-    }
-}
-
-impl<R: Real> Material<R> for SpirvMBFerromagnet<R> {
-    fn step_b(&mut self, context: &CellState<R>, delta_b: Vec3<R>) {
-        let target_b = context.with_m(self.m()).b() + delta_b;
-        let curve = coremem_types::mat::MBPgram::new(
-            self.0.curve.b_start,
-            self.0.curve.b_end,
-            self.0.curve.max_m,
-        );
-        // println!("step_b {}", target_b);
-        self.0.m = Vec3::new(
-            curve.move_b(self.0.m.x(), target_b.x()),
-            curve.move_b(self.0.m.y(), target_b.y()),
-            curve.move_b(self.0.m.z(), target_b.z()),
-        );
-    }
-    fn m(&self) -> Vec3<R> {
-        self.0.m()
-    }
-}
-
-// XXX: for debugging, use the same MBFerromagnet impl as we do in spirv impl.
-// pub type MBFerromagnet<R> = SpirvMBFerromagnet<R>;
-pub type MBFerromagnet<R> = NativeMBFerromagnet<R>;
-
-#[cfg(test)]
-mod test {
-    use super::*;
-    use float_eq::assert_float_eq;
-
-    #[test]
-    fn curve_interior() {
-        let curve = MBPgram::new(4.0, 6.0, 20.0f32);
-
-        assert_float_eq!(curve.move_b(0.0, 2.0), 0.0, abs <= 1e-5);
-        assert_float_eq!(curve.move_b(0.0, 5.0), 0.0, abs <= 1e-5);
-        assert_float_eq!(curve.move_b(1.0, 5.0), 1.0, abs <= 1e-5);
-        assert_float_eq!(curve.move_b(20.0, 5.0), 20.0, abs <= 1e-5);
-        assert_float_eq!(curve.move_b(-20.0, 4.0), -20.0, abs <= 1e-5);
-        assert_float_eq!(curve.move_b(-20.0, -6.0), -20.0, abs <= 1e-5);
-        assert_float_eq!(curve.move_b(20.0, -4.0), 20.0, abs <= 1e-5);
-        assert_float_eq!(curve.move_b(10.0, -2.0), 10.0, abs <= 1e-5);
-    }
-
-    #[test]
-    fn curve_exterior() {
-        let curve = MBPgram::new(4.0, 6.0, 20.0f32);
-        assert_float_eq!(curve.move_b(0.0, 6.0), 20.0, abs <= 1e-5);
-        assert_float_eq!(curve.move_b(0.0, 7.0), 20.0, abs <= 1e-5);
-        assert_float_eq!(curve.move_b(0.0, -6.0), -20.0, abs <= 1e-5);
-        assert_float_eq!(curve.move_b(0.0, -7.0), -20.0, abs <= 1e-5);
-
-
-        assert_float_eq!(curve.move_b(2.0, -6.0), -20.0, abs <= 1e-5);
-        assert_float_eq!(curve.move_b(20.0, -6.0), -20.0, abs <= 1e-5);
-        assert_float_eq!(curve.move_b(20.0, -5.0), 0.0, abs <= 1e-5);
-        assert_float_eq!(curve.move_b(20.0, -4.5), 10.0, abs <= 1e-5);
-        assert_float_eq!(curve.move_b(-15.0, 4.5), -10.0, abs <= 1e-5);
-        assert_float_eq!(curve.move_b(-15.0, 5.0), 0.0, abs <= 1e-5);
-        assert_float_eq!(curve.move_b(-15.0, 5.5), 10.0, abs <= 1e-5);
-        assert_float_eq!(curve.move_b(-15.0, 7.5), 20.0, abs <= 1e-5);
-    }
-
-    #[test]
-    fn curve_3r1() {
-        // slope of 3r1 is about M=793210*B
-        // This is almost identical to iron (795615)!
-        let curve = MBPgram::new(-0.3899, 0.3900, 310000f32);
-
-        // magnetizing:
-        // v.s. 198893 in B(H) curve
-        assert_float_eq!(curve.move_b(0.0, 0.250), 198703.0, abs <= 1.0);
-        // v.s. 278321 in B(H) curve
-        assert_float_eq!(curve.move_b(198703.0, 0.350), 278201.0, abs <= 1.0);
-        assert_float_eq!(curve.move_b(278201.0, 0.390), 310000.0, abs <= 1.0);
-
-        // de-magnetizing:
-        // From saturation, decreasing B causes NO decrease in M: instead, it causes a decrease in
-        // H. This is probably BAD: in the B(H) curve, a large change in H always causes a large
-        // change in B. The movement of H here is likely to induce current, whereas it SHOULDN'T.
-        assert_float_eq!(curve.move_b(310000.0, 0.38995), 310000.0, abs <= 1.0);
-        // v.s. 258626 in B(H); H = 220
-        assert_float_eq!(curve.move_b(310000.0, 0.325), 258406.0, abs <= 1.0);
-        // here's where H crosses 0 (v.s. B=0.325 in the B(H) curve... quite a difference)
-        assert_float_eq!(curve.move_b(310000.0, 0.050), 39788.438, abs <= 1.0);
-        // v.s. 35.0 in B(H)
-        assert_float_eq!(curve.move_b(310000.0, 0.0), 39.75, abs <= 1.0);
-
-        // negative magnetization:
-        // erase the magnetization: H = -40
-        assert_float_eq!(curve.move_b(310000.0, -0.00005), 0.0, abs <= 0.1);
-        // the magnetization has been completely erased:
-        assert_float_eq!(curve.move_b(310000.0, -0.25), -198703.0, abs <= 1.0);
-    }
-}

crates/coremem/src/meas.rs

@@ -1,80 +1,214 @@
-use crate::geom::{Meters, Region, Torus, Vec3, WorldRegion};
+use crate::geom::{HasCrossSection, Meters, Region, Torus, WorldRegion};
 use crate::real::{Real as _, ToFloat as _};
-use crate::sim::SampleableSim;
-use dyn_clone::{self, DynClone};
-use indexmap::IndexMap;
+use crate::cross::vec::{Vec3, Vec3u};
+use crate::sim::AbstractSim;
 use serde::{Serialize, Deserialize};
+use std::ops::AddAssign;
 
-// TODO: remove this Clone and Send requirement? Have Measurements be shared by-reference across
-// threads? i.e. Sync, and no Clone
-#[typetag::serde(tag = "type")]
-pub trait AbstractMeasurement: Send + Sync + DynClone {
-    fn eval(&self, state: &dyn SampleableSim) -> String;
-    fn key_value(&self, state: &dyn SampleableSim) -> IndexMap<String, String>;
+// TODO: do we really need both Send and Sync?
+pub trait AbstractMeasurement<S>: Send + Sync {
+    fn key_value(&self, state: &S) -> Vec<Measurement>;
 }
-dyn_clone::clone_trait_object!(AbstractMeasurement);
 
-pub fn eval_multiple_kv(state: &dyn SampleableSim, meas: &[Box<dyn AbstractMeasurement>]) -> IndexMap<String, String> {
-    let mut r = IndexMap::new();
-    for m in meas {
-        let other = m.key_value(state);
-        r.extend(other.into_iter());
+pub fn as_dyn_measurements<S, M: AbstractMeasurement<S>>(meas: &[M]) -> Vec<&dyn AbstractMeasurement<S>> {
+    meas.into_iter().map(|m| m as &dyn AbstractMeasurement<S>).collect()
+}
+
+
+/// combine several measurements
+pub fn eval_multiple<S>(state: &S, meas: &[&dyn AbstractMeasurement<S>]) -> Vec<Measurement> {
+    meas.into_iter().flat_map(|m| m.key_value(state).into_iter()).collect()
+}
+
+#[derive(Clone, Copy, Debug, PartialEq, Serialize, Deserialize)]
+pub enum MeasurementValue {
+    Field(Vec3<f32>),
+    Float(f32),
+    Int(u64),
+    Dim(Vec3u),
+}
+
+impl From<Vec3<f32>> for MeasurementValue {
+    fn from(v: Vec3<f32>) -> Self {
+        Self::Field(v)
+    }
+}
+impl From<f32> for MeasurementValue {
+    fn from(v: f32) -> Self {
+        Self::Float(v)
+    }
+}
+impl From<u64> for MeasurementValue {
+    fn from(v: u64) -> Self {
+        Self::Int(v)
+    }
+}
+impl From<Vec3u> for MeasurementValue {
+    fn from(v: Vec3u) -> Self {
+        Self::Dim(v)
+    }
+}
+
+#[derive(Clone, Debug, PartialEq, Serialize, Deserialize)]
+pub struct Measurement {
+    name: String,
+    value: MeasurementValue,
+    /// e.g. "A" for Amps
+    unit: String,
+}
+
+impl Measurement {
+    fn new<T: Into<MeasurementValue>>(name: &str, value: T, unit: &str) -> Self {
+        Self {
+            name: name.to_owned(),
+            value: value.into(),
+            unit: unit.to_owned(),
+        }
+    }
+    fn new_unitless<T: Into<MeasurementValue>>(name: &str, value: T) -> Self {
+        Self::new(name, value, "")
+    }
+
+    pub fn name(&self) -> &str {
+        &self.name
+    }
+
+    pub fn pretty_print(&self) -> String {
+        use MeasurementValue::*;
+        match self.value {
+            Field(v) => format!("{}{}", v, self.unit),
+            Float(f) => if self.unit != "" {
+                SiScale::format_short(f, &self.unit)
+            } else {
+                f.to_string()
+            },
+            Int(u)   => format!("{}{}", u, self.unit),
+            Dim(v)   => format!("{}x{}x{}{}", v.x(), v.y(), v.z(), self.unit),
+        }
+    }
+
+    /// format the Measurement in a way that could be parseable later.
+    /// one major use case for this is in dumping the type to a CSV.
+    pub fn machine_readable(&self) -> String {
+        use MeasurementValue::*;
+        match self.value {
+            Field(v) => format!("{},{},{}", v.x(), v.y(), v.z()),
+            Float(f) => f.to_string(),
+            Int(u)   => u.to_string(),
+            Dim(v)   => format!("{},{},{}", v.x(), v.y(), v.z()),
+        }
+    }
+
+    /// retrieve the float value of this measurement -- if it's of float type.
+    /// useful for tests
+    pub fn get_float(&self) -> Option<f32> {
+        match self.value {
+            MeasurementValue::Float(f) => Some(f),
+            _ => None,
+        }
+    }
+}
+
+impl<S> AbstractMeasurement<S> for Measurement {
+    fn key_value(&self, _state: &S) -> Vec<Measurement> {
+        vec![self.clone()]
+    }
+}
+
+enum SiScale {
+    Pico,
+    Nano,
+    Micro,
+    Milli,
+    Unit,
+    Kilo,
+    Mega,
+    Giga,
+    Terra,
+}
+
+impl SiScale {
+    fn for_value(v: f32) -> Self {
+        use SiScale::*;
+        match v.abs() {
+            v if v < 1e-12 => Unit,
+            v if v < 1e-9  => Pico,
+            v if v < 1e-6  => Nano,
+            v if v < 1e-3  => Micro,
+            v if v < 1e0   => Milli,
+            v if v < 1e3   => Unit,
+            v if v < 1e6   => Kilo,
+            v if v < 1e9   => Mega,
+            v if v < 1e12  => Giga,
+            v if v < 1e15  => Terra,
+            _ => Unit
+        }
+    }
+
+    /// return the numerical scale of this prefix.
+    /// e.g. `scale(&Pico) -> 1e-12
+    fn scale(&self) -> f32 {
+        use SiScale::*;
+        match *self {
+            Pico =>  1e-12,
+            Nano =>  1e-9,
+            Micro => 1e-6,
+            Milli => 1e-3,
+            Unit =>  1.0,
+            Kilo =>  1e3,
+            Mega =>  1e6,
+            Giga =>  1e9,
+            Terra => 1e12,
+        }
+    }
+
+    /// return the short string for this scale.
+    /// e.g. `shortcode(Pico) -> "p"`
+    fn shortcode(&self) -> &'static str {
+        use SiScale::*;
+        match *self {
+            Pico =>  "p",
+            Nano =>  "n",
+            Micro => "u",
+            Milli => "m",
+            Unit =>  "",
+            Kilo =>  "k",
+            Mega =>  "M",
+            Giga =>  "G",
+            Terra => "T",
+        }
+    }
+
+    /// format `v`, with the provided unit.
+    /// e.g. `format_short(1234, "A") -> "1.23 kA"
+    fn format_short(v: f32, unit: &str) -> String {
+        let si = SiScale::for_value(v);
+        let scaled = v / si.scale();
+        format!("{:.2} {}{}", scaled, si.shortcode(), unit)
     }
-    r
 }
 
 #[derive(Clone, Serialize, Deserialize)]
 pub struct Time;
 
-#[typetag::serde]
-impl AbstractMeasurement for Time {
-    fn eval(&self, state: &dyn SampleableSim) -> String {
-        format!("{:.3e}s (step {})", state.time(), state.step_no())
-    }
-    fn key_value(&self, state: &dyn SampleableSim) -> IndexMap<String, String> {
-        [
-            ("step".to_string(), state.step_no().to_string()),
-            ("time".to_string(), state.time().to_string()),
-        ].into_iter().collect()
+impl<S: AbstractSim> AbstractMeasurement<S> for Time {
+    fn key_value(&self, state: &S) -> Vec<Measurement> {
+        vec![
+            Measurement::new_unitless("step", state.step_no()),
+            Measurement::new("time", state.time(), "s"),
+        ]
     }
 }
 
 #[derive(Clone, Serialize, Deserialize)]
 pub struct Meta;
 
-#[typetag::serde]
-impl AbstractMeasurement for Meta {
-    fn eval(&self, state: &dyn SampleableSim) -> String {
-        format!("{}x{}x{} feat: {:.1e}m", state.width(), state.height(), state.depth(), state.feature_size())
-    }
-    fn key_value(&self, state: &dyn SampleableSim) -> IndexMap<String, String> {
-        [
-            ("width".to_string(), state.width().to_string()),
-            ("height".to_string(), state.height().to_string()),
-            ("depth".to_string(), state.depth().to_string()),
-            ("feature_size".to_string(), state.feature_size().to_string()),
-        ].into_iter().collect()
-    }
-}
-
-#[derive(Clone, Serialize, Deserialize)]
-pub struct Label(pub String);
-
-impl Label {
-    pub fn new<S: Into<String>>(s: S) -> Self {
-        Self(s.into())
-    }
-}
-
-#[typetag::serde]
-impl AbstractMeasurement for Label {
-    fn eval(&self, _state: &dyn SampleableSim) -> String {
-        self.0.clone()
-    }
-    fn key_value(&self, _state: &dyn SampleableSim) -> IndexMap<String, String> {
-        [
-            (self.0.clone(), self.0.clone()),
-        ].into_iter().collect()
+impl<S: AbstractSim> AbstractMeasurement<S> for Meta {
+    fn key_value(&self, state: &S) -> Vec<Measurement> {
+        vec![
+            Measurement::new_unitless("dim", state.size().0),
+            Measurement::new("feature_size", state.feature_size(), "m"),
+        ]
     }
 }
 
@@ -92,25 +226,18 @@ impl Volume {
         }
     }
     /// Returns the volume of the region, in units of um^3
-    fn data(&self, state: &dyn SampleableSim) -> f32 {
+    fn data<S: AbstractSim>(&self, state: &S) -> f32 {
         let feat_um = state.feature_size() as f64 * 1e6;
         
         (state.volume_of_region(&*self.region) as f64 * feat_um * feat_um * feat_um) as f32
     }
 }
 
-#[typetag::serde]
-impl AbstractMeasurement for Volume {
-    fn eval(&self, state: &dyn SampleableSim) -> String {
-        format!("Vol({}): {:.2e} um^3",
-            self.name,
-            self.data(state),
-        )
-    }
-    fn key_value(&self, state: &dyn SampleableSim) -> IndexMap<String, String> {
-        [
-            (format!("Vol({})", self.name), self.data(state).to_string()),
-        ].into_iter().collect()
+impl<S: AbstractSim> AbstractMeasurement<S> for Volume {
+    fn key_value(&self, state: &S) -> Vec<Measurement> {
+        vec![
+            Measurement::new(&format!("Vol({})", self.name), self.data(state), "um^3"),
+        ]
     }
 }
 
@@ -127,7 +254,7 @@ impl Current {
             region: Box::new(r)
         }
     }
-    fn data(&self, state: &dyn SampleableSim) -> (f32, Vec3<f32>) {
+    fn data<S: AbstractSim>(&self, state: &S) -> (f32, Vec3<f32>) {
         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().cast(), current.cast())
@@ -138,6 +265,24 @@ impl Current {
     }
 }
 
+// TODO: clean up these FieldSample types
+#[derive(Default)]
+struct TupleSum<T>(T);
+
+impl<T0: Default + AddAssign, T1: Default + AddAssign> std::iter::Sum for TupleSum<(T0, T1)> {
+    fn sum<I>(iter: I) -> Self
+        where I: Iterator<Item = Self>
+    {
+        let mut s = Self::default();
+        for TupleSum((a0, a1)) in iter {
+            s.0.0 += a0;
+            s.0.1 += a1;
+        }
+        s
+    }
+}
+
+
 #[derive(Default)]
 struct FieldSample(u32, f64, Vec3<f64>);
 
@@ -183,65 +328,79 @@ impl std::iter::Sum for FieldSamples<[FieldSample; 3]> {
     }
 }
 
-#[typetag::serde]
-impl AbstractMeasurement for Current {
-    fn eval(&self, state: &dyn SampleableSim) -> String {
+impl<S: AbstractSim> AbstractMeasurement<S> for Current {
+    fn key_value(&self, state: &S) -> Vec<Measurement> {
         let (mean_current_mag, mean_current_vec) = self.data(state);
-        format!("I/cell({}): {:.2e} {:.2e}",
-            self.name,
-            mean_current_mag,
-            mean_current_vec)
-    }
-    fn key_value(&self, state: &dyn SampleableSim) -> IndexMap<String, String> {
-        let (mean_current_mag, mean_current_vec) = self.data(state);
-        [
-            (format!("Imag/cell({})", self.name), mean_current_mag.to_string()),
-            (format!("I/cell({})", self.name), mean_current_vec.to_string()),
-        ].into_iter().collect()
+        vec![
+            Measurement::new(
+                &format!("Imag/cell({})", self.name),
+                mean_current_mag,
+                "A",
+            ),
+            Measurement::new(
+                &format!("I/cell({})", self.name),
+                mean_current_vec,
+                "A",
+            ),
+        ]
     }
 }
 
 /// Measures the current directed around a closed loop
 #[derive(Clone, Serialize, Deserialize)]
-pub struct CurrentLoop {
+pub struct CurrentLoop<R> {
     name: String,
-    region: Torus
+    region: R,
 }
 
-impl CurrentLoop {
-    pub fn new(name: &str, r: Torus) -> Self {
+impl<R> CurrentLoop<R> {
+    pub fn new(name: &str, r: R) -> Self {
         Self {
             name: name.into(),
             region: r,
         }
     }
-    fn data(&self, state: &dyn SampleableSim) -> 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.cast());
-            FieldSample(1, directed_current.cast(), current.cast())
+}
+impl<R: Region + HasCrossSection> CurrentLoop<R> {
+    fn data<S: AbstractSim>(&self, state: &S) -> f32 {
+        // - current exists as a property of a 2d surface.
+        // - the user has provided us a 3d volume which behaves as though it's an extruded surface:
+        //   for any point in the volume we can query the normal vector of the cross section
+        //   containing that point.
+        // - we choose that measuring the "current" on such a volume means to measure the average
+        //   current through all its cross sections (for most boring materials, each
+        //   cross section has nearly identical current).
+        // - therefore, enumerate the entire volume and compute the "net" current (the sum over
+        //   each cell of whatever current in that cell is along the cross-section normal).
+        //   then divide by the number of complete cross sections we measured, to average.
+        let feature_area = state.feature_size() * state.feature_size();
+        let TupleSum((net_current, cross_sections)) = state.map_sum_over_enumerated(&self.region, move |coord: Meters, _cell| {
+            // `normal` represents both the size of the cross section (m^2) this cell belongs to,
+            // and the normal direction of the cross section.
+            let normal = self.region.cross_section_normal(coord);  // [m^2]
+            // calculate the amount of normal current through this specific cell
+            let current_density = state.current_density(coord);    // [A/m^2]
+            let cross_sectional_current = feature_area * current_density.dot(normal.norm()); // [A]
+            // keep track of how many cross sections we enumerate, since each additional cross
+            // sections represents a double-count of the current.
+            let num_cross_sections_filled = feature_area / normal.mag();
+            TupleSum((cross_sectional_current, num_cross_sections_filled))
         });
-        let mean_directed_current = directed_current.cast::<f32>() / 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
+        let mean_cross_sectional_current = net_current.cast::<f32>() / cross_sections;
+        mean_cross_sectional_current
     }
 }
 
-#[typetag::serde]
-impl AbstractMeasurement for CurrentLoop {
-    fn eval(&self, state: &dyn SampleableSim) -> String {
+impl<R: Region + HasCrossSection, S: AbstractSim> AbstractMeasurement<S> for CurrentLoop<R> {
+    fn key_value(&self, state: &S) -> Vec<Measurement> {
         let cross_sectional_current = self.data(state);
-        format!("I({}): {:.2e}", self.name, cross_sectional_current)
-    }
-    fn key_value(&self, state: &dyn SampleableSim) -> IndexMap<String, String> {
-        let cross_sectional_current = self.data(state);
-        [
-            (format!("I({})", self.name), cross_sectional_current.to_string()),
-        ].into_iter().collect()
+        vec![
+            Measurement::new(
+                &format!("I({})", self.name),
+                cross_sectional_current,
+                "A"
+            ),
+        ]
     }
 }
 
@@ -260,7 +419,7 @@ impl MagneticLoop {
             region: r,
         }
     }
-    fn data(&self, state: &dyn SampleableSim) -> (f32, f32, f32) {
+    fn data<S: AbstractSim>(&self, state: &S) -> (f32, f32, f32) {
         let FieldSamples([
             FieldSample(volume, directed_m, _m_vec),
             FieldSample(_, directed_b, _b_vec),
@@ -300,24 +459,14 @@ impl MagneticLoop {
     }
 }
 
-#[typetag::serde]
-impl AbstractMeasurement for MagneticLoop {
-    fn eval(&self, state: &dyn SampleableSim) -> String {
+impl<S: AbstractSim> AbstractMeasurement<S> for MagneticLoop {
+    fn key_value(&self, state: &S) -> Vec<Measurement> {
         let (mean_directed_m, mean_directed_b, mean_directed_h) = self.data(state);
-        format!(
-            "M({}): {:.2e}; B({}): {:.2e}; H({}): {:.2e}",
-            self.name, mean_directed_m,
-            self.name, mean_directed_b,
-            self.name, mean_directed_h,
-        )
-    }
-    fn key_value(&self, state: &dyn SampleableSim) -> IndexMap<String, String> {
-        let (mean_directed_m, mean_directed_b, mean_directed_h) = self.data(state);
-        [
-            (format!("M({})", self.name), mean_directed_m.to_string()),
-            (format!("B({})", self.name), mean_directed_b.to_string()),
-            (format!("H({})", self.name), mean_directed_h.to_string()),
-        ].into_iter().collect()
+        vec![
+            Measurement::new_unitless(&format!("M({})", self.name), mean_directed_m),
+            Measurement::new_unitless(&format!("B({})", self.name), mean_directed_b),
+            Measurement::new_unitless(&format!("H({})", self.name), mean_directed_h),
+        ]
     }
 }
 
@@ -335,7 +484,7 @@ impl MagneticFlux {
             region: Box::new(r)
         }
     }
-    fn data(&self, state: &dyn SampleableSim) -> Vec3<f32> {
+    fn data<S: AbstractSim>(&self, state: &S) -> Vec3<f32> {
         let FieldSample(volume, _directed_mag, mag_vec) = state.map_sum_over(&*self.region, |cell| {
             let b = cell.b();
             let mag = b.mag();
@@ -346,17 +495,15 @@ impl MagneticFlux {
     }
 }
 
-#[typetag::serde]
-impl AbstractMeasurement for MagneticFlux {
-    fn eval(&self, state: &dyn SampleableSim) -> String {
+impl<S: AbstractSim> AbstractMeasurement<S> for MagneticFlux {
+    fn key_value(&self, state: &S) -> Vec<Measurement> {
         let mean_mag = self.data(state);
-        format!("Bavg({}): {:.2e}", self.name, mean_mag)
-    }
-    fn key_value(&self, state: &dyn SampleableSim) -> IndexMap<String, String> {
-        let mean_mag = self.data(state);
-        [
-            (format!("Bavg({})", self.name), mean_mag.to_string()),
-        ].into_iter().collect()
+        vec![
+            Measurement::new_unitless(
+                &format!("Bavg({})", self.name),
+                mean_mag,
+            )
+        ]
     }
 }
 
@@ -374,7 +521,7 @@ impl Magnetization {
             region: Box::new(r)
         }
     }
-    fn data(&self, state: &dyn SampleableSim) -> Vec3<f32> {
+    fn data<S: AbstractSim>(&self, state: &S) -> Vec3<f32> {
         let FieldSample(volume, _directed_mag, mag_vec) = state.map_sum_over(&*self.region, |cell| {
             let m = cell.m();
             let mag = m.mag();
@@ -385,17 +532,14 @@ impl Magnetization {
     }
 }
 
-#[typetag::serde]
-impl AbstractMeasurement for Magnetization {
-    fn eval(&self, state: &dyn SampleableSim) -> String {
+impl<S: AbstractSim> AbstractMeasurement<S> for Magnetization {
+    fn key_value(&self, state: &S) -> Vec<Measurement> {
         let mean_mag = self.data(state);
-        format!("Mavg({}): {:.2e}", self.name, mean_mag)
-    }
-    fn key_value(&self, state: &dyn SampleableSim) -> IndexMap<String, String> {
-        let mean_mag = self.data(state);
-        [
-            (format!("Mavg({})", self.name), mean_mag.to_string()),
-        ].into_iter().collect()
+        vec![
+            Measurement::new_unitless(
+                &format!("Mavg({})", self.name), mean_mag
+            ),
+        ]
     }
 }
 
@@ -407,17 +551,12 @@ fn loc(v: Meters) -> String {
 #[derive(Clone, Serialize, Deserialize)]
 pub struct MagnetizationAt(pub Meters);
 
-#[typetag::serde]
-impl AbstractMeasurement for MagnetizationAt {
-    fn eval(&self, state: &dyn SampleableSim) -> String {
+impl<S: AbstractSim> AbstractMeasurement<S> for MagnetizationAt {
+    fn key_value(&self, state: &S) -> Vec<Measurement> {
         let m = state.sample(self.0).m();
-        format!("M{}: {:.2e}", loc(self.0), m)
-    }
-    fn key_value(&self, state: &dyn SampleableSim) -> IndexMap<String, String> {
-        let m = state.sample(self.0).m();
-        [
-            (format!("M{}", loc(self.0)), m.to_string()),
-        ].into_iter().collect()
+        vec![
+            Measurement::new_unitless(&format!("M{}", loc(self.0)), m.cast())
+        ]
     }
 }
 
@@ -425,17 +564,14 @@ impl AbstractMeasurement for MagnetizationAt {
 #[derive(Clone, Serialize, Deserialize)]
 pub struct MagneticFluxAt(pub Meters);
 
-#[typetag::serde]
-impl AbstractMeasurement for MagneticFluxAt {
-    fn eval(&self, state: &dyn SampleableSim) -> String {
+impl<S: AbstractSim> AbstractMeasurement<S> for MagneticFluxAt {
+    fn key_value(&self, state: &S) -> Vec<Measurement> {
         let b = state.sample(self.0).b();
-        format!("B{}: {:.2e}", loc(self.0), b)
-    }
-    fn key_value(&self, state: &dyn SampleableSim) -> IndexMap<String, String> {
-        let b = state.sample(self.0).b();
-        [
-            (format!("B{}", loc(self.0)), b.to_string()),
-        ].into_iter().collect()
+        vec![
+            Measurement::new_unitless(
+                &format!("B{}", loc(self.0)), b.cast()
+            )
+        ]
     }
 }
 
@@ -443,34 +579,28 @@ impl AbstractMeasurement for MagneticFluxAt {
 #[derive(Clone, Serialize, Deserialize)]
 pub struct MagneticStrengthAt(pub Meters);
 
-#[typetag::serde]
-impl AbstractMeasurement for MagneticStrengthAt {
-    fn eval(&self, state: &dyn SampleableSim) -> String {
+impl<S: AbstractSim> AbstractMeasurement<S> for MagneticStrengthAt {
+    fn key_value(&self, state: &S) -> Vec<Measurement> {
         let h = state.sample(self.0).h();
-        format!("H{}: {:.2e}", loc(self.0), h)
-    }
-    fn key_value(&self, state: &dyn SampleableSim) -> IndexMap<String, String> {
-        let h = state.sample(self.0).h();
-        [
-            (format!("H{}", loc(self.0)), h.to_string()),
-        ].into_iter().collect()
+        vec![
+            Measurement::new_unitless(
+                &format!("H{}", loc(self.0)), h.cast()
+            )
+        ]
     }
 }
 
 #[derive(Clone, Serialize, Deserialize)]
 pub struct ElectricField(pub Meters);
 
-#[typetag::serde]
-impl AbstractMeasurement for ElectricField {
-    fn eval(&self, state: &dyn SampleableSim) -> String {
+impl<S: AbstractSim> AbstractMeasurement<S> for ElectricField {
+    fn key_value(&self, state: &S) -> Vec<Measurement> {
         let e = state.sample(self.0).e();
-        format!("E{}: {:.2e}", loc(self.0), e)
-    }
-    fn key_value(&self, state: &dyn SampleableSim) -> IndexMap<String, String> {
-        let e = state.sample(self.0).e();
-        [
-            (format!("E{}", loc(self.0)), e.to_string()),
-        ].into_iter().collect()
+        vec![
+            Measurement::new_unitless(
+                &format!("E{}", loc(self.0)), e.cast()
+            )
+        ]
     }
 }
 
@@ -490,7 +620,7 @@ impl Energy {
             region: Box::new(region),
         }
     }
-    fn data(&self, state: &dyn SampleableSim) -> f32 {
+    fn data<S: AbstractSim>(&self, state: &S) -> 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
@@ -507,17 +637,14 @@ impl Energy {
     }
 }
 
-#[typetag::serde]
-impl AbstractMeasurement for Energy {
-    fn eval(&self, state: &dyn SampleableSim) -> String {
+impl<S: AbstractSim> AbstractMeasurement<S> for Energy {
+    fn key_value(&self, state: &S) -> Vec<Measurement> {
         let e = self.data(state);
-        format!("U({}): {:.2e}", self.name, e)
-    }
-    fn key_value(&self, state: &dyn SampleableSim) -> IndexMap<String, String> {
-        let e = self.data(state);
-        [
-            (format!("U({})", self.name), e.to_string()),
-        ].into_iter().collect()
+        vec![
+            Measurement::new(
+                &format!("U({})", self.name), e, "J"
+            )
+        ]
     }
 }
 
@@ -537,7 +664,7 @@ impl Power {
             region: Box::new(region),
         }
     }
-    fn data(&self, state: &dyn SampleableSim) -> f32 {
+    fn data<S: AbstractSim>(&self, state: &S) -> 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)]
@@ -549,16 +676,191 @@ impl Power {
     }
 }
 
-#[typetag::serde]
-impl AbstractMeasurement for Power {
-    fn eval(&self, state: &dyn SampleableSim) -> String {
+impl<S: AbstractSim> AbstractMeasurement<S> for Power {
+    fn key_value(&self, state: &S) -> Vec<Measurement> {
         let power = self.data(state);
-        format!("P({}): {:.2e}", self.name, power)
-    }
-    fn key_value(&self, state: &dyn SampleableSim) -> IndexMap<String, String> {
-        let power = self.data(state);
-        [
-            (format!("P({})", self.name), power.to_string()),
-        ].into_iter().collect()
+        vec![
+            Measurement::new(
+                &format!("P({})", self.name), power, "W"
+            )
+        ]
+    }
+}
+
+#[cfg(test)]
+pub mod test {
+    use super::*;
+    use crate::cross::mat::AnisomorphicConductor;
+    use crate::cross::step::SimMeta;
+    use crate::geom::Index;
+    use crate::sim::{Fields, GenericSim, StaticSim};
+    use crate::stim::AbstractStimulus;
+
+    struct MockSim {
+        e_field: Vec3<f32>,
+        dim: Vec3u,
+        feature_size: f32,
+        mat: AnisomorphicConductor<f32>,
+    }
+    impl AbstractSim for MockSim {
+        type Real = f32;
+        type Material = AnisomorphicConductor<f32>;
+
+        fn meta(&self) -> SimMeta<f32> {
+            SimMeta::new(self.dim, self.feature_size, 1e-9)
+        }
+        fn step_no(&self) -> u64 {
+            unimplemented!()
+        }
+        fn fields_at_index(&self, _pos: Index) -> Fields<Self::Real> {
+            Fields::new(self.e_field, Vec3::zero(), Vec3::zero())
+        }
+        fn get_material_index(&self, _at: Index) -> &Self::Material {
+            &self.mat
+        }
+        fn put_material_index(&mut self, _at: Index, _m: Self::Material) {
+            unimplemented!()
+        }
+        fn step_multiple<S: AbstractStimulus>(&mut self, _num_steps: u32, _s: &S) {
+            unimplemented!()
+        }
+        fn to_static(&self) -> StaticSim {
+            unimplemented!()
+        }
+        fn to_generic(&self) -> GenericSim<Self::Real> {
+            unimplemented!()
+        }
+    }
+
+    #[derive(Clone, Serialize, Deserialize)]
+    struct MockRegion {
+        normal: Vec3<f32>,
+    }
+    impl HasCrossSection for MockRegion {
+        fn cross_section_normal(&self, _p: Meters) -> Vec3<f32> {
+            self.normal
+        }
+    }
+    #[typetag::serde]
+    impl Region for MockRegion {
+        fn contains(&self, _p: Meters) -> bool {
+            true
+        }
+    }
+
+    #[test]
+    fn current_loop_trivial() {
+        let sim = MockSim {
+            e_field: Vec3::new(1.0, 0.0, 0.0),
+            dim: Vec3u::new(1, 1, 1),
+            feature_size: 1.0,
+            mat: AnisomorphicConductor::new(Vec3::new(1.0, 1.0, 1.0)),
+        };
+        let region = MockRegion {
+            normal: Vec3::new(1.0, 0.0, 0.0),
+        };
+
+        let kv = CurrentLoop::new("test", region).key_value(&sim);
+        assert_eq!(kv.len(), 1);
+        // measured area is 1 m^2
+        // region cross-section is 1 m^2
+        // conductivity is 1 S/m
+        assert_eq!(kv[0].get_float().unwrap(), 1.0);
+    }
+    #[test]
+    fn current_loop_multi_cell() {
+        let sim = MockSim {
+            e_field: Vec3::new(1.0, 0.0, 0.0),
+            dim: Vec3u::new(4, 4, 4),
+            feature_size: 0.25,
+            mat: AnisomorphicConductor::new(Vec3::new(1.0, 1.0, 1.0)),
+        };
+        let region = MockRegion {
+            normal: Vec3::new(1.0, 0.0, 0.0),
+        };
+
+        let kv = CurrentLoop::new("test", region).key_value(&sim);
+        assert_eq!(kv.len(), 1);
+        // measured area is 1 m^2
+        // region cross-section is 1 m^2
+        // conductivity is 1 S/m
+        assert_eq!(kv[0].get_float().unwrap(), 1.0);
+    }
+    #[test]
+    fn current_loop_off_conductor() {
+        let sim = MockSim {
+            e_field: Vec3::new(1.0, 1.0, 1.0),
+            dim: Vec3u::new(4, 4, 4),
+            feature_size: 0.25,
+            mat: AnisomorphicConductor::new(Vec3::new(0.0, 1.0, 1.0)),
+        };
+        let region = MockRegion {
+            normal: Vec3::new(1.0, 0.0, 0.0),
+        };
+
+        let kv = CurrentLoop::new("test", region).key_value(&sim);
+        assert_eq!(kv.len(), 1);
+        // material is not conductive in the direction being queried
+        assert_eq!(kv[0].get_float().unwrap(), 0.0);
+    }
+    #[test]
+    fn current_loop_e_field() {
+        let sim = MockSim {
+            e_field: Vec3::new(4.0, 2.0, 1.0),
+            dim: Vec3u::new(4, 4, 4),
+            feature_size: 0.25,
+            mat: AnisomorphicConductor::new(Vec3::new(1.0, 1.0, 1.0)),
+        };
+        let region = MockRegion {
+            normal: Vec3::new(1.0, 0.0, 0.0),
+        };
+
+        let kv = CurrentLoop::new("test", region).key_value(&sim);
+        assert_eq!(kv.len(), 1);
+        // measured area is 1 m^2
+        // region cross-section is 1 m^2
+        // conductivity is 1 S/m
+        // e field is 4 V/m
+        assert_eq!(kv[0].get_float().unwrap(), 4.0);
+    }
+    #[test]
+    fn current_loop_conductivity() {
+        let sim = MockSim {
+            e_field: Vec3::new(4.0, 2.0, 1.0),
+            dim: Vec3u::new(4, 4, 4),
+            feature_size: 0.25,
+            mat: AnisomorphicConductor::new(Vec3::new(3.0, 1.0, 1.0)),
+        };
+        let region = MockRegion {
+            normal: Vec3::new(1.0, 0.0, 0.0),
+        };
+
+        let kv = CurrentLoop::new("test", region).key_value(&sim);
+        assert_eq!(kv.len(), 1);
+        // measured area is 1 m^2
+        // region cross-section is 1 m^2
+        // conductivity is 3 S/m
+        // e field is 4 V/m
+        assert_eq!(kv[0].get_float().unwrap(), 3.0*4.0);
+    }
+    #[test]
+    fn current_loop_cross_section() {
+        let sim = MockSim {
+            e_field: Vec3::new(4.0, 2.0, 1.0),
+            dim: Vec3u::new(4, 4, 4),
+            feature_size: 0.5,
+            mat: AnisomorphicConductor::new(Vec3::new(3.0, 1.0, 1.0)),
+        };
+        let region = MockRegion {
+            normal: Vec3::new(16.0, 0.0, 0.0),
+        };
+
+        let kv = CurrentLoop::new("test", region).key_value(&sim);
+        assert_eq!(kv.len(), 1);
+        // measured area is 2 m^2
+        // region cross-section is 16 m^2
+        // conductivity is 3 S/m
+        // e field is 4 V/m
+        assert_eq!(kv[0].get_float().unwrap(), 3.0*4.0*16.0);
     }
 }

crates/coremem/src/render.rs

@@ -1,7 +1,8 @@
-use crate::geom::{Meters, Vec2, Vec3};
+use crate::geom::Index;
 use crate::real::ToFloat as _;
-use crate::sim::{SampleableSim, Sample, StaticSim};
-use crate::meas::{self, AbstractMeasurement};
+use crate::cross::vec::{Vec2, Vec3};
+use crate::sim::{AbstractSim, GenericSim, Sample};
+use crate::meas::{self, AbstractMeasurement, Measurement};
 use crossterm::{cursor, QueueableCommand as _};
 use crossterm::style::{style, Color, PrintStyledContent, Stylize as _};
 use font8x8::{BASIC_FONTS, GREEK_FONTS, UnicodeFonts as _};
@@ -11,6 +12,8 @@ use image::{RgbImage, Rgb};
 use imageproc::{pixelops, drawing};
 use rayon::prelude::*;
 use serde::{Serialize, Deserialize};
+use std::collections::hash_map::DefaultHasher;
+use std::hash::Hasher;
 use std::fs::{File, OpenOptions};
 use std::io::{BufReader, BufWriter, Seek as _, SeekFrom, Write as _};
 use std::path::{Path, PathBuf};
@@ -51,10 +54,10 @@ fn scale_unsigned_to_u8(x: f32, typ: f32) -> u8 {
 /// Scale a vector to have magnitude between [0, 1).
 fn scale_vector(x: Vec2<f32>, typical_mag: f32) -> Vec2<f32> {
     let new_mag = scale_unsigned(x.mag(), typical_mag);
-    x.with_mag(new_mag)
+    x.with_mag(new_mag).unwrap_or_default()
 }
 
-fn im_size<S: SampleableSim>(state: &S, max_w: u32, max_h: u32) -> (u32, u32) {
+fn im_size<S: AbstractSim>(state: &S, max_w: u32, max_h: u32) -> (u32, u32) {
     let mut width = max_w;
     let mut height = width * state.height() / state.width();
     if height > max_h {
@@ -71,6 +74,7 @@ pub enum FieldDisplayMode {
     EzBxy,
     BCurrent,
     M,
+    Material,
 }
 
 impl FieldDisplayMode {
@@ -80,17 +84,19 @@ impl FieldDisplayMode {
             BzExy => EzBxy,
             EzBxy => BCurrent,
             BCurrent => M,
-            M => BzExy,
+            M => Material,
+            Material => BzExy,
         }
     }
 
     pub fn prev(self) -> Self {
         use FieldDisplayMode::*;
         match self {
-            BzExy => M,
+            BzExy => Material,
             EzBxy => BzExy,
             BCurrent => EzBxy,
             M => BCurrent,
+            Material => M,
         }
     }
 }
@@ -128,20 +134,22 @@ impl RenderConfig {
 struct RenderSteps<'a, S> {
     im: RgbImage, 
     sim: &'a S,
-    meas: &'a [Box<dyn AbstractMeasurement>],
+    meas: &'a [&'a dyn AbstractMeasurement<S>],
     /// Simulation z coordinate to sample
     z: u32,
 }
 
-impl<'a, S: SampleableSim> RenderSteps<'a, S> {
+impl<'a, S: AbstractSim> RenderSteps<'a, S> {
+    // TODO: this could probably be a single measurement, and we just let collections of
+    // measurements also behave as measurements
     /// Render using default configuration constants
-    fn render(state: &'a S, measurements: &'a [Box<dyn AbstractMeasurement>], z: u32) -> RgbImage {
+    fn render(state: &'a S, measurements: &'a [&'a dyn AbstractMeasurement<S>], z: u32) -> RgbImage {
         Self::render_configured(state, measurements, z, (640, 480), RenderConfig::default())
     }
     /// Render, controlling things like the size.
     fn render_configured(
         state: &'a S,
-        measurements: &'a [Box<dyn AbstractMeasurement>],
+        measurements: &'a [&'a dyn AbstractMeasurement<S>],
         z: u32,
         max_size: (u32, u32),
         config: RenderConfig,
@@ -173,11 +181,14 @@ impl<'a, S: SampleableSim> RenderSteps<'a, S> {
             FieldDisplayMode::M => {
                 me.render_m(config.scale);
             }
+            FieldDisplayMode::Material => {
+                me.render_mat(config.scale);
+            }
         }
         me.render_measurements();
         me.im
     }
-    fn new(sim: &'a S, meas: &'a [Box<dyn AbstractMeasurement>], width: u32, height: u32, z: u32) -> Self {
+    fn new(sim: &'a S, meas: &'a [&'a dyn AbstractMeasurement<S>], width: u32, height: u32, z: u32) -> Self {
         RenderSteps {
             im: RgbImage::new(width, height),
             sim,
@@ -186,13 +197,10 @@ impl<'a, S: SampleableSim> RenderSteps<'a, S> {
         }
     }
 
-    fn get_at_px(&self, x_px: u32, y_px: u32) -> Sample {
-        let x_prop = x_px as f32 / self.im.width() as f32;
-        let x_m = x_prop * (self.sim.width() as f32 * self.sim.feature_size() as f32);
-        let y_prop = y_px as f32 / self.im.height() as f32;
-        let y_m = y_prop * (self.sim.height() as f32 * self.sim.feature_size() as f32);
-        let z_m = self.z as f32 * self.sim.feature_size() as f32;
-        self.sim.sample(Meters(Vec3::new(x_m, y_m, z_m)))
+    fn get_at_px<'b>(&'b self, x_px: u32, y_px: u32) -> Sample<'b, S::Real, S::Material> {
+        let x_idx = x_px * self.sim.width() / self.im.width();
+        let y_idx = y_px * self.sim.height() / self.im.height();
+        self.sim.sample(Index::new(x_idx, y_idx, self.z))
     }
 
     //////////////   Ex/Ey/Bz configuration ////////////
@@ -229,7 +237,22 @@ impl<'a, S: SampleableSim> RenderSteps<'a, S> {
         self.render_vector_field(Rgb([0xff, 0xff, 0xff]), 1.0e5 * scale, |cell| cell.m().xy().to_f32());
     }
 
-    fn render_vector_field<F: Fn(&Sample) -> Vec2<f32>>(&mut self, color: Rgb<u8>, typical: f32, measure: F) {
+    fn render_mat(&mut self, scale: f32) {
+        unsafe fn to_bytes<T>(d: &T) -> &[u8] {
+            std::slice::from_raw_parts(d as *const T as *const u8, std::mem::size_of::<T>())
+        }
+        self.render_scalar_field(scale, false, 1, |cell| {
+            let mut hasher = DefaultHasher::new();
+            let as_bytes = unsafe { to_bytes(cell.material()) };
+            std::hash::Hash::hash_slice(as_bytes, &mut hasher);
+            hasher.finish() as f32 / (-1i64 as u64 as f32)
+        });
+    }
+
+    fn render_vector_field<F>(&mut self, color: Rgb<u8>, typical: f32, measure: F)
+    where
+        F: Fn(&Sample<'_, S::Real, S::Material>) -> Vec2<f32>
+    {
         let w = self.im.width();
         let h = self.im.height();
         let vec_spacing = 10;
@@ -244,7 +267,10 @@ impl<'a, S: SampleableSim> RenderSteps<'a, S> {
             }
         }
     }
-    fn render_scalar_field<F: Fn(&Sample) -> f32 + Sync>(&mut self, typical: f32, signed: bool, slot: u32, measure: F) {
+    fn render_scalar_field<F>(&mut self, typical: f32, signed: bool, slot: u32, measure: F)
+    where
+        F: Fn(&Sample<'_, S::Real, S::Material>) -> f32 + Sync
+    {
         // XXX: get_at_px borrows self, so we need to clone the image to operate on it mutably.
         let mut im = self.im.clone();
         let w = im.width();
@@ -268,8 +294,8 @@ impl<'a, S: SampleableSim> RenderSteps<'a, S> {
         self.im = im;
     }
     fn render_measurements(&mut self) {
-        for (meas_no, m) in self.meas.iter().enumerate() {
-            let meas_string = m.eval(self.sim);
+        for (meas_no, m) in meas::eval_multiple(self.sim, &self.meas).into_iter().enumerate() {
+            let meas_string = m.pretty_print();
             for (i, c) in meas_string.chars().enumerate() {
                 let glyph = BASIC_FONTS.get(c)
                     .or_else(|| GREEK_FONTS.get(c))
@@ -290,7 +316,10 @@ impl<'a, S: SampleableSim> RenderSteps<'a, S> {
         }
     }
 
-    fn field_vector<F: Fn(&Sample) -> Vec2<f32>>(&self, xidx: u32, yidx: u32, size: u32, measure: &F) -> Vec2<f32> {
+    fn field_vector<F>(&self, xidx: u32, yidx: u32, size: u32, measure: &F) -> Vec2<f32>
+    where
+        F: Fn(&Sample<'_, S::Real, S::Material>) -> Vec2<f32>
+    {
         let mut field = Vec2::default();
         let w = self.im.width();
         let h = self.im.height();
@@ -339,25 +368,25 @@ impl ImageRenderExt for RgbImage {
 }
 
 pub trait Renderer<S>: Send + Sync {
-    fn render_z_slice(&self, state: &S, z: u32, measurements: &[Box<dyn AbstractMeasurement>], config: RenderConfig);
+    fn render_z_slice(&self, state: &S, z: u32, measurements: &[&dyn AbstractMeasurement<S>], config: RenderConfig);
     // {
     //     self.render_with_image(state, &RenderSteps::render(state, measurements, z), measurements);
     // }
-    fn render(&self, state: &S, measurements: &[Box<dyn AbstractMeasurement>], config: RenderConfig);
+    fn render(&self, state: &S, measurements: &[&dyn AbstractMeasurement<S>], config: RenderConfig);
     /// Not intended to be called directly by users; implement this if you want the image to be
     /// computed using default settings and you just manage where to display/save it.
-    fn render_with_image(&self, state: &S, _im: &RgbImage, measurements: &[Box<dyn AbstractMeasurement>], config: RenderConfig) {
+    fn render_with_image(&self, state: &S, _im: &RgbImage, measurements: &[&dyn AbstractMeasurement<S>], config: RenderConfig) {
         self.render(state, measurements, config);
     }
 }
 
-fn default_render_z_slice<S: SampleableSim, R: Renderer<S>>(
-    me: &R, state: &S, z: u32, measurements: &[Box<dyn AbstractMeasurement>], config: RenderConfig,
+fn default_render_z_slice<S: AbstractSim, R: Renderer<S>>(
+    me: &R, state: &S, z: u32, measurements: &[&dyn AbstractMeasurement<S>], config: RenderConfig,
 ) {
     me.render_with_image(state, &RenderSteps::render(state, measurements, z), measurements, config);
 }
-fn default_render<S: SampleableSim, R: Renderer<S>>(
-    me: &R, state: &S, measurements: &[Box<dyn AbstractMeasurement>], config: RenderConfig
+fn default_render<S: AbstractSim, R: Renderer<S>>(
+    me: &R, state: &S, measurements: &[&dyn AbstractMeasurement<S>], config: RenderConfig
 ) {
         me.render_z_slice(state, state.depth() / 2, measurements, config);
 }
@@ -365,7 +394,7 @@ fn default_render<S: SampleableSim, R: Renderer<S>>(
 // pub struct NumericTermRenderer;
 // 
 // impl Renderer for NumericTermRenderer {
-//     fn render(&mut self, state: &SimSnapshot, _measurements: &[Box<dyn AbstractMeasurement>]) {
+//     fn render(&mut self, state: &SimSnapshot, _measurements: &[&dyn AbstractMeasurement<S>]) {
 //         for y in 0..state.height() {
 //             for x in 0..state.width() {
 //                 let cell = state.get((x, y).into());
@@ -385,17 +414,18 @@ fn default_render<S: SampleableSim, R: Renderer<S>>(
 #[derive(Default)]
 pub struct ColorTermRenderer;
 
-impl<S: SampleableSim> Renderer<S> for ColorTermRenderer {
-    fn render(&self, state: &S, measurements: &[Box<dyn AbstractMeasurement>], config: RenderConfig) {
+impl<S: AbstractSim> Renderer<S> for ColorTermRenderer {
+    fn render(&self, state: &S, measurements: &[&dyn AbstractMeasurement<S>], config: RenderConfig) {
         default_render(self, state, measurements, config)
     }
     fn render_z_slice(
         &self,
         state: &S,
         z: u32,
-        measurements: &[Box<dyn AbstractMeasurement>],
+        measurements: &[&dyn AbstractMeasurement<S>],
         config: RenderConfig,
     ) {
+        let measurements = meas::eval_multiple(state, measurements);
         let (max_w, mut max_h) = crossterm::terminal::size().unwrap();
         max_h = max_h.saturating_sub(2 + measurements.len() as u16);
         let im = RenderSteps::render_configured(state, &[], z, (max_w as _, max_h as _), config);
@@ -424,7 +454,7 @@ impl<S: SampleableSim> Renderer<S> for ColorTermRenderer {
         for m in measurements {
             // Measurements can be slow to compute
             stdout.flush().unwrap();
-            let meas_string = m.eval(state);
+            let meas_string = m.pretty_print();
             stdout.queue(cursor::MoveDown(1)).unwrap();
             stdout.queue(cursor::MoveToColumn(1)).unwrap();
             stdout.queue(PrintStyledContent(style(meas_string))).unwrap();
@@ -447,14 +477,14 @@ impl Y4MRenderer {
     }
 }
 
-impl<S: SampleableSim> Renderer<S> for Y4MRenderer {
-    fn render_z_slice(&self, state: &S, z: u32, measurements: &[Box<dyn AbstractMeasurement>], config: RenderConfig) {
+impl<S: AbstractSim> Renderer<S> for Y4MRenderer {
+    fn render_z_slice(&self, state: &S, z: u32, measurements: &[&dyn AbstractMeasurement<S>], config: RenderConfig) {
         default_render_z_slice(self, state, z, measurements, config)
     }
-    fn render(&self, state: &S, measurements: &[Box<dyn AbstractMeasurement>], config: RenderConfig) {
+    fn render(&self, state: &S, measurements: &[&dyn AbstractMeasurement<S>], config: RenderConfig) {
         default_render(self, state, measurements, config)
     }
-    fn render_with_image(&self, _state: &S, im: &RgbImage, _meas: &[Box<dyn AbstractMeasurement>], _config: RenderConfig) {
+    fn render_with_image(&self, _state: &S, im: &RgbImage, _meas: &[&dyn AbstractMeasurement<S>], _config: RenderConfig) {
         {
             let mut enc = self.encoder.lock().unwrap();
             if enc.is_none() {
@@ -539,17 +569,17 @@ impl<S> MultiRenderer<S> {
     }
 }
 
-impl<S: SampleableSim> Renderer<S> for MultiRenderer<S> {
-    fn render_z_slice(&self, state: &S, z: u32, measurements: &[Box<dyn AbstractMeasurement>], config: RenderConfig) {
+impl<S: AbstractSim> Renderer<S> for MultiRenderer<S> {
+    fn render_z_slice(&self, state: &S, z: u32, measurements: &[&dyn AbstractMeasurement<S>], config: RenderConfig) {
         default_render_z_slice(self, state, z, measurements, config)
     }
-    fn render(&self, state: &S, measurements: &[Box<dyn AbstractMeasurement>], config: RenderConfig) {
+    fn render(&self, state: &S, measurements: &[&dyn AbstractMeasurement<S>], config: RenderConfig) {
         if self.renderers.read().unwrap().len() != 0 {
             self.render_with_image(state, &RenderSteps::render(state, measurements, state.depth() / 2), measurements, config);
         }
     }
 
-    fn render_with_image(&self, state: &S, im: &RgbImage, measurements: &[Box<dyn AbstractMeasurement>], config: RenderConfig) {
+    fn render_with_image(&self, state: &S, im: &RgbImage, measurements: &[&dyn AbstractMeasurement<S>], config: RenderConfig) {
         for r in &*self.renderers.read().unwrap() {
             if r.work_this_frame(state.step_no()) {
                 r.renderer.render_with_image(state, im, measurements, config);
@@ -559,25 +589,28 @@ impl<S: SampleableSim> Renderer<S> for MultiRenderer<S> {
 }
 
 #[derive(Serialize, Deserialize)]
-pub struct SerializedFrame<S=StaticSim> {
+pub struct SerializedFrame<S> {
     pub state: S,
     /// although not generally necessary to load the sim, saving the measurements is beneficial for
     /// post-processing.
-    pub measurements: Vec<Box<dyn AbstractMeasurement>>,
+    pub measurements: Vec<Measurement>,
 }
 
-impl<S: SampleableSim> SerializedFrame<S> {
-    pub fn to_static(self) -> SerializedFrame<StaticSim> {
+impl<S: AbstractSim> SerializedFrame<S> {
+    pub fn to_generic(self) -> SerializedFrame<GenericSim<S::Real>> {
         SerializedFrame {
-            state: SampleableSim::to_static(&self.state),
+            state: AbstractSim::to_generic(&self.state),
             measurements: self.measurements,
         }
     }
 }
 
+/// this serializes the simulation state plus measurements to disk.
+/// it can either convert the state to a generic, material-agnostic format (generic)
+/// or dump it as-is.
 pub struct SerializerRenderer {
     fmt_str: String,
-    prefer_static: bool,
+    prefer_generic: bool,
 }
 
 impl SerializerRenderer {
@@ -586,25 +619,25 @@ impl SerializerRenderer {
     pub fn new(fmt_str: &str) -> Self {
         Self {
             fmt_str: fmt_str.into(),
-            prefer_static: false,
+            prefer_generic: false,
         }
     }
 
-    /// Same as `new`, but cast to StaticSim before serializing. This yields a file that's easier
-    /// for post-processing, and may be smaller in size.
-    pub fn new_static(fmt_str: &str) -> Self {
+    /// Same as `new`, but cast to GenericSim before serializing. This yields a file that's easier
+    /// for post-processing.
+    pub fn new_generic(fmt_str: &str) -> Self {
         Self {
             fmt_str: fmt_str.into(),
-            prefer_static: true,
+            prefer_generic: true,
         }
     }
 }
 
 impl SerializerRenderer {
-    fn serialize<S: SampleableSim + Serialize>(&self, state: &S, measurements: &[Box<dyn AbstractMeasurement>]) {
+    fn serialize<S: AbstractSim + Serialize>(&self, state: &S, measurements: Vec<Measurement>) {
         let frame = SerializedFrame {
             state,
-            measurements: measurements.iter().cloned().collect(),
+            measurements,
         };
         let name = self.fmt_str.replace("{step_no}", &*frame.state.step_no().to_string());
         let out = BufWriter::new(File::create(name).unwrap());
@@ -612,21 +645,21 @@ impl SerializerRenderer {
         bincode::serialize_into(out, &frame).unwrap();
     }
 
-    pub fn try_load<S: SampleableSim + for <'a> Deserialize<'a>>(&self) -> Option<SerializedFrame<S>> {
+    pub fn try_load<S: AbstractSim + for <'a> Deserialize<'a>>(&self) -> Option<SerializedFrame<S>> {
         let mut reader = BufReader::new(File::open(&*self.fmt_str).ok()?);
         bincode::deserialize_from(&mut reader).ok()
     }
 }
 
-impl<S: SampleableSim + Serialize> Renderer<S> for SerializerRenderer {
-    fn render_z_slice(&self, state: &S, z: u32, measurements: &[Box<dyn AbstractMeasurement>], config: RenderConfig) {
+impl<S: AbstractSim + Serialize> Renderer<S> for SerializerRenderer {
+    fn render_z_slice(&self, state: &S, z: u32, measurements: &[&dyn AbstractMeasurement<S>], config: RenderConfig) {
         default_render_z_slice(self, state, z, measurements, config)
     }
-    fn render(&self, state: &S, measurements: &[Box<dyn AbstractMeasurement>], _config: RenderConfig) {
-        if self.prefer_static {
-            self.serialize(&state.to_static(), measurements);
+    fn render(&self, state: &S, measurements: &[&dyn AbstractMeasurement<S>], _config: RenderConfig) {
+        if self.prefer_generic {
+            self.serialize(&state.to_generic(), meas::eval_multiple(state, measurements));
         } else {
-            self.serialize(state, measurements);
+            self.serialize(state, meas::eval_multiple(state, measurements));
         }
     }
 }
@@ -661,12 +694,12 @@ impl CsvRenderer {
     }
 }
 
-impl<S: SampleableSim> Renderer<S> for CsvRenderer {
-    fn render_z_slice(&self, state: &S, z: u32, measurements: &[Box<dyn AbstractMeasurement>], config: RenderConfig) {
+impl<S: AbstractSim> Renderer<S> for CsvRenderer {
+    fn render_z_slice(&self, state: &S, z: u32, measurements: &[&dyn AbstractMeasurement<S>], config: RenderConfig) {
         default_render_z_slice(self, state, z, measurements, config)
     }
-    fn render(&self, state: &S, measurements: &[Box<dyn AbstractMeasurement>], _config: RenderConfig) {
-        let row = meas::eval_multiple_kv(state, measurements);
+    fn render(&self, state: &S, measurements: &[&dyn AbstractMeasurement<S>], _config: RenderConfig) {
+        let row = meas::eval_multiple(state, measurements);
         let step = state.step_no();
         let mut lock = self.state.lock().unwrap();
         let mut writer = match lock.take().unwrap() {
@@ -700,13 +733,13 @@ impl<S: SampleableSim> Renderer<S> for CsvRenderer {
                     file.set_len(0).unwrap();
                     let mut writer = csv::Writer::from_writer(BufWriter::new(file));
                     // write the header
-                    writer.write_record(row.keys()).unwrap();
+                    writer.write_record(row.iter().map(|m| m.name())).unwrap();
                     writer
                 }
             },
             CsvState::Writing(writer) => writer,
         };
-        writer.write_record(row.values()).unwrap();
+        writer.write_record(row.iter().map(|m| m.machine_readable())).unwrap();
         writer.flush().unwrap();
         *lock = Some(CsvState::Writing(writer));
     }

crates/coremem/src/sim/legacy/mat/bh_ferromagnet.rs

@@ -1,8 +1,9 @@
-use crate::CellState;
-use crate::geom::{Line2d, Vec2, Vec3, Polygon2d};
-use crate::mat::Material;
+use super::Material;
+use crate::material_compat;
+use crate::geom::{Line2d, Polygon2d};
 use crate::real::Real;
-use crate::sim::StepParametersMut;
+use crate::sim::legacy::{CellState, StepParametersMut};
+use crate::cross::vec::{Vec2, Vec3};
 
 use lazy_static::lazy_static;
 use log::trace;
@@ -199,6 +200,7 @@ impl<R: Real> Material<R> for Ferroxcube3R1<R> {
         StepParametersMut::default().with_conductivity(Vec3::uniform(1e-3))
     }
 }
+material_compat!(R, Ferroxcube3R1<R>);
 
 /// Simple, square-loop ferrite
 #[derive(Default, Copy, Clone, PartialEq, Serialize, Deserialize)]
@@ -235,6 +237,7 @@ impl<R: Real> Material<R> for MinimalSquare<R> {
         StepParametersMut::default().with_conductivity(Vec3::uniform(1e-3))
     }
 }
+material_compat!(R, MinimalSquare<R>);
 
 
 #[cfg(test)]

crates/coremem/src/sim/legacy/mat/db.rs

@@ -1,8 +1,8 @@
 //! database of common materials
 
-use crate::geom::Vec3;
-use crate::mat::{AnisomorphicConductor, IsomorphicConductor, LinearMagnet, Ferroxcube3R1, MinimalSquare};
+use super::{AnisomorphicConductor, IsomorphicConductor, LinearMagnet, Ferroxcube3R1, MinimalSquare};
 use crate::real::Real;
+use crate::cross::vec::Vec3;
 
 pub fn conductor<R: Real, R2: Real>(conductivity: R2) -> IsomorphicConductor<R> {
     IsomorphicConductor::new(conductivity.cast())

crates/coremem/src/sim/legacy/mat/linear.rs

@@ -1,7 +1,8 @@
-use crate::CellState;
-use crate::geom::Vec3;
-use crate::mat::Material;
+use super::Material;
+use crate::material_compat;
 use crate::real::Real;
+use crate::sim::legacy::CellState;
+use crate::cross::vec::Vec3;
 
 use serde::{Serialize, Deserialize};
 
@@ -45,6 +46,7 @@ impl<R: Real> Material<R> for LinearMagnet<R> {
         self.m += delta_m;
     }
 }
+material_compat!(R, LinearMagnet<R>);
 
 #[cfg(test)]
 mod test {

crates/coremem/src/sim/legacy/mat/mod.rs

@@ -1,29 +1,24 @@
-use crate::CellState;
-use crate::geom::Vec3;
 use crate::real::Real;
-use crate::sim::{PmlParameters, PmlState, StepParameters, StepParametersMut};
+use crate::sim::legacy::{CellState, PmlParameters, PmlState, StepParameters, StepParametersMut};
+use crate::cross::vec::Vec3;
 
-use enum_dispatch::enum_dispatch;
 use serde::{Serialize, Deserialize};
 
 pub mod db;
 mod bh_ferromagnet;
-mod mb_ferromagnet;
 mod linear;
 
 pub use bh_ferromagnet::*;
-pub use mb_ferromagnet::*;
-pub use coremem_types::mat::{
+pub use coremem_cross::mat::{
     AnisomorphicConductor,
     Ferroxcube3R1MH,
     FullyGenericMaterial,
     IsoConductorOr,
     IsomorphicConductor,
-    MHPgram
+    MHPgram,
 };
 pub use linear::*;
 
-#[enum_dispatch]
 pub trait Material<R: Real> {
     fn step_parameters_mut<'a>(&'a mut self) -> StepParametersMut<'a, R> {
         // by default, behave as a vacuum
@@ -37,6 +32,22 @@ pub trait Material<R: Real> {
     fn step_b(&mut self, _context: &CellState<R>, _delta_b: Vec3<R>) {
     }
 }
+#[macro_export]
+macro_rules! material_compat {
+    (R, $mat:path) => {
+        // XXX this is not that useful an implementation.
+        // it exists mostly because some users want the `Material::conductivity()` method.
+        impl<R: Real> crate::mat::Material<R> for $mat {
+            fn conductivity(&self) -> Vec3<R> {
+                crate::sim::legacy::mat::MaterialExt::step_parameters(self).conductivity()
+            }
+            fn move_b_vec(&self, _m: Vec3<R>, _target_b: Vec3<R>) -> Vec3<R> {
+                unimplemented!()
+            }
+        }
+    }
+}
+
 
 pub trait MaterialExt<R> {
     fn step_parameters<'a>(&'a self) -> StepParameters<'a, R>;
@@ -87,6 +98,7 @@ impl<R: Real> Material<R> for Static<R> {
         self.m
     }
 }
+material_compat!(R, Static<R>);
 
 impl<R: Real, T> From<T> for Static<R>
 where T: Into<GenericMaterial<R>>
@@ -111,15 +123,15 @@ impl<R: Real> Material<R> for Pml<R> {
         StepParametersMut::default().with_pml(&mut self.0, self.1)
     }
 }
+material_compat!(R, Pml<R>);
 
 
-// #[enum_dispatch(Material)]
 #[derive(Clone, PartialEq, Serialize, Deserialize)]
 pub enum GenericMaterial<R> {
     Conductor(AnisomorphicConductor<R>),
     LinearMagnet(LinearMagnet<R>),
     Pml(Pml<R>),
-    MBFerromagnet(MBFerromagnet<R>),
+    MBPgram(MBPgram<R>),
     Ferroxcube3R1(Ferroxcube3R1<R>),
     MinimalSquare(MinimalSquare<R>),
 }
@@ -155,9 +167,9 @@ impl<R> From<Pml<R>> for GenericMaterial<R> {
     }
 }
 
-impl<R> From<MBFerromagnet<R>> for GenericMaterial<R> {
-    fn from(inner: MBFerromagnet<R>) -> Self {
-        Self::MBFerromagnet(inner)
+impl<R> From<MBPgram<R>> for GenericMaterial<R> {
+    fn from(inner: MBPgram<R>) -> Self {
+        Self::MBPgram(inner)
     }
 }
 
@@ -180,7 +192,7 @@ impl<R: Real> Material<R> for GenericMaterial<R> {
             Conductor(inner) => inner.step_parameters_mut(),
             LinearMagnet(inner) => inner.step_parameters_mut(),
             Pml(inner) => inner.step_parameters_mut(),
-            MBFerromagnet(inner) => inner.step_parameters_mut(),
+            MBPgram(inner) => inner.step_parameters_mut(),
             Ferroxcube3R1(inner) => inner.step_parameters_mut(),
             MinimalSquare(inner) => inner.step_parameters_mut(),
         }
@@ -192,7 +204,7 @@ impl<R: Real> Material<R> for GenericMaterial<R> {
             Conductor(inner) => inner.m(),
             LinearMagnet(inner) => inner.m(),
             Pml(inner) => inner.m(),
-            MBFerromagnet(inner) => inner.m(),
+            MBPgram(inner) => inner.m(),
             Ferroxcube3R1(inner) => Material::m(inner),
             MinimalSquare(inner) => Material::m(inner),
         }
@@ -204,19 +216,19 @@ impl<R: Real> Material<R> for GenericMaterial<R> {
             Conductor(inner) => inner.step_b(context, delta_b),
             LinearMagnet(inner) => inner.step_b(context, delta_b),
             Pml(inner) => inner.step_b(context, delta_b),
-            MBFerromagnet(inner) => inner.step_b(context, delta_b),
+            MBPgram(inner) => inner.step_b(context, delta_b),
             Ferroxcube3R1(inner) => inner.step_b(context, delta_b),
             MinimalSquare(inner) => inner.step_b(context, delta_b),
         }
     }
 }
+material_compat!(R, GenericMaterial<R>);
 
-// #[enum_dispatch(Material)]
 #[derive(Clone, Serialize, Deserialize)]
 pub enum GenericMaterialNoPml<R> {
     Conductor(AnisomorphicConductor<R>),
     LinearMagnet(LinearMagnet<R>),
-    MBFerromagnet(MBFerromagnet<R>),
+    MBPgram(MBPgram<R>),
     Ferroxcube3R1(Ferroxcube3R1<R>),
     MinimalSquare(MinimalSquare<R>),
 }
@@ -239,7 +251,7 @@ impl<R: Real> Material<R> for GenericMaterialNoPml<R> {
         match self {
             Conductor(inner) => inner.step_parameters_mut(),
             LinearMagnet(inner) => inner.step_parameters_mut(),
-            MBFerromagnet(inner) => inner.step_parameters_mut(),
+            MBPgram(inner) => inner.step_parameters_mut(),
             Ferroxcube3R1(inner) => inner.step_parameters_mut(),
             MinimalSquare(inner) => inner.step_parameters_mut(),
         }
@@ -250,7 +262,7 @@ impl<R: Real> Material<R> for GenericMaterialNoPml<R> {
         match self {
             Conductor(inner) => inner.m(),
             LinearMagnet(inner) => inner.m(),
-            MBFerromagnet(inner) => inner.m(),
+            MBPgram(inner) => inner.m(),
             Ferroxcube3R1(inner) => Material::m(inner),
             MinimalSquare(inner) => Material::m(inner),
         }
@@ -261,16 +273,16 @@ impl<R: Real> Material<R> for GenericMaterialNoPml<R> {
         match self {
             Conductor(inner) => inner.step_b(context, delta_b),
             LinearMagnet(inner) => inner.step_b(context, delta_b),
-            MBFerromagnet(inner) => inner.step_b(context, delta_b),
+            MBPgram(inner) => inner.step_b(context, delta_b),
             Ferroxcube3R1(inner) => inner.step_b(context, delta_b),
             MinimalSquare(inner) => inner.step_b(context, delta_b),
         }
     }
 }
+material_compat!(R, GenericMaterialNoPml<R>);
 
 
 /// Materials which have only 1 Vec3.
-// #[enum_dispatch(Material)]
 #[derive(Clone, Serialize, Deserialize)]
 pub enum GenericMaterialOneField<R> {
     Conductor(AnisomorphicConductor<R>),
@@ -318,17 +330,66 @@ impl<R: Real> Material<R> for GenericMaterialOneField<R> {
         }
     }
 }
+material_compat!(R, GenericMaterialOneField<R>);
 
-// coremem_types adapters
+// coremem_cross adapters
+// TODO: move this to a dedicated file
+
+/// the coremem_cross Materials are stateless;
+/// rather than hold onto their own magnetic fields (for example), the simulation holds that.
+/// that's counter to the original cpu-only simulation, in which materials hold their own state.
+///
+/// this type adapts any stateless coremem_cross::mat::Material type to be a coremem::mat::Material.
+#[derive(Default, Copy, Clone, PartialEq, Serialize, Deserialize)]
+pub struct AdaptStateless<R, M> {
+    mat: M,
+    m: Vec3<R>,
+}
+
+impl<R, M> AdaptStateless<R, M> {
+    pub fn into_inner(self) -> M {
+        self.mat
+    }
+}
+
+impl<R: Default, M> From<M> for AdaptStateless<R, M> {
+    fn from(mat: M) -> Self {
+        Self { mat, m: Default::default() }
+    }
+}
+
+impl<R: Real, M: coremem_cross::mat::Material<R>> Material<R> for AdaptStateless<R, M> {
+    fn m(&self) -> Vec3<R> {
+        self.m
+    }
+    fn step_parameters_mut<'a>(&'a mut self) -> StepParametersMut<'a, R> {
+        let c = self.mat.conductivity();
+        StepParametersMut::default().with_conductivity(c)
+    }
+    fn step_b(&mut self, context: &CellState<R>, delta_b: Vec3<R>) {
+        let target_b = context.with_m(self.m).b() + delta_b;
+        self.m = self.mat.move_b_vec(self.m, target_b);
+    }
+}
+
+// conductors: these are truly stateless
 impl<R: Real> Material<R> for AnisomorphicConductor<R> {
     fn step_parameters_mut<'a>(&'a mut self) -> StepParametersMut<'a, R> {
-        let c = coremem_types::mat::Material::conductivity(self);
+        let c = coremem_cross::mat::Material::conductivity(self);
         StepParametersMut::default().with_conductivity(c)
     }
 }
 impl<R: Real> Material<R> for IsomorphicConductor<R> {
     fn step_parameters_mut<'a>(&'a mut self) -> StepParametersMut<'a, R> {
-        let c = coremem_types::mat::Material::conductivity(self);
+        let c = coremem_cross::mat::Material::conductivity(self);
         StepParametersMut::default().with_conductivity(c)
     }
 }
+
+pub type MBPgram<R> = AdaptStateless<R, coremem_cross::mat::MBPgram<R>>;
+impl<R: Real> MBPgram<R> {
+    pub fn new(b_start: R, b_end: R, m_max: R) -> Self {
+        coremem_cross::mat::MBPgram::new(b_start, b_end, m_max).into()
+    }
+}
+material_compat!(R, MBPgram<R>);

crates/coremem/src/sim/spirv/bindings.rs

@@ -1,194 +0,0 @@
-use serde::de::Deserializer;
-use serde::ser::Serializer;
-use serde::{Deserialize, Serialize};
-
-use crate::mat::{
-    AnisomorphicConductor,
-    IsoConductorOr,
-    IsomorphicConductor,
-    Ferroxcube3R1MH,
-    FullyGenericMaterial,
-    MaterialExt as _,
-    MBFerromagnet,
-    MBPgram,
-    MHPgram,
-    Static
-};
-use crate::geom::{Index, Vec3, Vec3u};
-
-/// hide the actual spirv backend structures inside a submodule to make their use/boundary clear.
-mod ffi {
-    pub use spirv_backend::entry_points;
-    pub use spirv_backend::sim::SerializedSimMeta;
-    pub use spirv_backend::support::Optional;
-    pub use coremem_types::mat::MBPgram;
-}
-
-// conversion traits for types defined cross-lib
-pub trait IntoFfi {
-    type Ffi;
-    fn into_ffi(self) -> Self::Ffi;
-}
-
-pub trait IntoLib {
-    type Lib;
-    fn into_lib(self) -> Self::Lib;
-}
-
-macro_rules! identity {
-    ($($param:ident,)* => $t:ty) => {
-        impl<$($param: IntoFfi),*> IntoFfi for $t {
-            type Ffi = $t;
-            fn into_ffi(self) -> Self::Ffi {
-                self
-            }
-        }
-        impl<$($param: IntoLib),*> IntoLib for $t {
-            type Lib = $t;
-            fn into_lib(self) -> Self::Lib {
-                self
-            }
-        }
-    };
-}
-
-// XXX: should work for any other lifetime, not just 'static
-identity!(=> f32);
-identity!(=> &'static str);
-identity!(T0, T1, => (T0, T1));
-identity!(=> Vec3u);
-identity!(T, => Vec3<T>);
-
-impl<L: IntoFfi> IntoFfi for Option<L>
-    where L::Ffi: Default
-{
-    type Ffi = ffi::Optional<L::Ffi>;
-    fn into_ffi(self) -> Self::Ffi {
-        match self {
-            Some(s) => ffi::Optional::some(s.into_ffi()),
-            None => ffi::Optional::none(),
-        }
-    }
-}
-
-impl<F: Copy + IntoLib> IntoLib for ffi::Optional<F> {
-    type Lib = Option<F::Lib>;
-    fn into_lib(self) -> Self::Lib {
-        if self.is_some() {
-            Some(self.unwrap().into_lib())
-        } else {
-            None
-        }
-    }
-}
-
-impl IntoFfi for MBPgram<f32> {
-    type Ffi = ffi::MBPgram<f32>;
-    fn into_ffi(self) -> Self::Ffi {
-        Self::Ffi::new(self.b_start, self.b_end, self.max_m)
-    }
-}
-
-impl IntoLib for ffi::MBPgram<f32> {
-    type Lib = MBPgram<f32>;
-    fn into_lib(self) -> Self::Lib {
-        Self::Lib::new(self.b_start, self.b_end, self.max_m)
-    }
-}
-
-identity!( => MHPgram<f32>);
-identity!( => Ferroxcube3R1MH);
-identity!( => FullyGenericMaterial<f32>);
-identity!(R, M, => IsoConductorOr<R, M>);
-
-// N.B.: this isn't fully correct, as Static also encodes the magnetic component
-impl From<Static<f32>> for FullyGenericMaterial<f32> {
-    fn from(m: Static<f32>) -> Self {
-        AnisomorphicConductor::new(m.conductivity).into()
-    }
-}
-
-impl From<MBFerromagnet<f32>> for FullyGenericMaterial<f32> {
-    fn from(m: MBFerromagnet<f32>) -> Self {
-        m.curve().into_ffi().into()
-    }
-}
-
-// this is bitwise- and type-compatible with the spirv SimMeta, except we need serde traits
-#[derive(Clone, Default, Serialize, Deserialize)]
-pub struct SimMeta {
-    pub(crate) dim: Index,
-    pub(crate) inv_feature_size: f32,
-    pub(crate) time_step: f32,
-    pub(crate) feature_size: f32,
-}
-
-impl IntoFfi for SimMeta {
-    type Ffi = ffi::SerializedSimMeta;
-    fn into_ffi(self) -> Self::Ffi {
-        Self::Ffi {
-            dim: self.dim.0.into_ffi(),
-            inv_feature_size: self.inv_feature_size,
-            time_step: self.time_step,
-            feature_size: self.feature_size,
-        }
-    }
-}
-
-/// Store the FFI form in memory, but serialize via the lib form.
-#[derive(Clone, Default, PartialEq)]
-pub struct Remote<F>(F);
-
-impl<F> Remote<F> {
-    pub fn into_inner(self) -> F {
-        self.0
-    }
-}
-
-impl<L: IntoFfi> From<L> for Remote<L::Ffi> {
-    fn from(l: L) -> Self {
-        Remote(l.into_ffi())
-    }
-}
-
-impl<F> std::ops::Deref for Remote<F> {
-    type Target = F;
-    fn deref(&self) -> &Self::Target {
-        &self.0
-    }
-}
-
-impl<F> Serialize for Remote<F>
-where F: Clone + IntoLib,
-      F::Lib: Serialize,
-{
-    fn serialize<S>(&self, serializer: S) -> Result<S::Ok, S::Error>
-    where
-        S: Serializer,
-    {
-        let local = self.0.clone().into_lib();
-        local.serialize(serializer)
-    }
-}
-
-impl<'de, F> Deserialize<'de> for Remote<F>
-    where F: IntoLib,
-          F::Lib: Deserialize<'de> + IntoFfi<Ffi=F>,
-{
-    fn deserialize<D>(deserializer: D) -> Result<Self, D::Error>
-    where
-        D: Deserializer<'de>,
-    {
-        let local: F::Lib = Deserialize::deserialize(deserializer)?;
-        Ok(Remote(local.into_ffi()))
-    }
-}
-
-// FUNCTION BINDINGS
-pub fn entry_points<L>() -> Option<(&'static str, &'static str)>
-where
-    L: IntoFfi,
-    L::Ffi: 'static
-{
-    ffi::entry_points::<L::Ffi>().into_lib()
-}

crates/coremem/src/sim/spirv/cpu.rs

@@ -0,0 +1,46 @@
+use coremem_cross::mat::Material;
+use coremem_cross::real::Real;
+use coremem_cross::step::{SimMeta, StepEContext, StepHContext};
+use coremem_cross::vec::{Vec3, Vec3u};
+
+use super::SimBackend;
+
+#[derive(Default)]
+pub struct CpuBackend;
+
+impl<R: Real, M: Material<R>> SimBackend<R, M> for CpuBackend {
+    fn step_n(
+        &mut self,
+        meta: SimMeta<R>,
+        mat: &[M],
+        stim_e: &[Vec3<R>],
+        stim_h: &[Vec3<R>],
+        e: &mut [Vec3<R>],
+        h: &mut [Vec3<R>],
+        m: &mut [Vec3<R>],
+        num_steps: u32,
+    ) {
+        for _ in 0..num_steps {
+            // step E field
+            apply_all_cells(meta.dim, |idx| {
+                StepEContext::step_flat_view(meta, mat, stim_e, e, h, idx);
+            });
+
+            // step H field
+            apply_all_cells(meta.dim, |idx| {
+                StepHContext::step_flat_view(meta, mat, stim_h, e, h, m, idx);
+            });
+        }
+    }
+}
+
+fn apply_all_cells<F: FnMut(Vec3u)>(dim: Vec3u, mut f: F) {
+    for z in 0..dim.z() {
+        for y in 0..dim.y() {
+            for x in 0..dim.x() {
+                f(Vec3u::new(x, y, z));
+            }
+        }
+    }
+}
+

crates/coremem/src/sim/spirv/gpu.rs

@@ -0,0 +1,444 @@
+use futures::FutureExt as _;
+use log::info;
+use std::borrow::Cow;
+use std::num::NonZeroU64;
+use wgpu;
+use wgpu::util::DeviceExt as _;
+
+use coremem_cross::vec::{Vec3, Vec3u};
+use coremem_cross::step::SimMeta;
+
+use spirv_backend::HasEntryPoints;
+
+use super::SimBackend;
+
+#[derive(Default)]
+pub struct WgpuBackend {
+    handles: Option<(&'static str /* step_h */, &'static str /* step_e */, WgpuHandles)>,
+}
+
+struct WgpuHandles {
+    step_bind_group_layout: wgpu::BindGroupLayout,
+    step_e_pipeline: wgpu::ComputePipeline,
+    step_h_pipeline: wgpu::ComputePipeline,
+    device: wgpu::Device,
+    queue: wgpu::Queue,
+}
+
+impl WgpuHandles {
+    fn open<R, M: HasEntryPoints<R>>(dim: Vec3u) -> Self {
+        info!("WgpuHandles::open({})", dim);
+        use std::mem::size_of;
+        let volume = dim.product_sum_usize() as u64;
+        let max_elem_size = size_of::<M>().max(size_of::<Vec3<R>>());
+        let max_array_size = volume * max_elem_size as u64;
+        let max_buf_size = max_array_size + 0x1000; // allow some overhead
+
+        let (device, queue) = futures::executor::block_on(open_device(max_buf_size));
+        let shader_binary = get_shader();
+        let shader_module = unsafe { device.create_shader_module_spirv(&shader_binary) };
+        let (step_bind_group_layout, step_h_pipeline, step_e_pipeline) = make_pipelines(
+            &device, &shader_module, M::step_h(), M::step_e()
+        );
+        WgpuHandles {
+            step_bind_group_layout,
+            step_h_pipeline,
+            step_e_pipeline,
+            device,
+            queue,
+        }
+    }
+}
+
+// TODO: these bounds aren't 100% right. we're sending R and M over to the GPU by a bitwise copy.
+// that probably means the types should be Send + Copy
+impl<R: Copy, M: Send + Sync + HasEntryPoints<R>> SimBackend<R, M> for WgpuBackend {
+    fn step_n(
+        &mut self,
+        meta: SimMeta<R>,
+        mat: &[M],
+        stim_cpu_e: &[Vec3<R>],
+        stim_cpu_h: &[Vec3<R>],
+        e: &mut [Vec3<R>],
+        h: &mut [Vec3<R>],
+        m: &mut [Vec3<R>],
+        num_steps: u32,
+    ) {
+        let field_bytes = meta.dim.product_sum() as usize * std::mem::size_of::<Vec3<f32>>();
+
+        let (step_h, step_e, handles) = self.handles.get_or_insert_with(|| (
+            M::step_h(),
+            M::step_e(),
+            WgpuHandles::open::<R, M>(meta.dim)
+        ));
+        // if device is opened, make sure we're open for the right types
+        assert_eq!(*step_h, M::step_h());
+        assert_eq!(*step_e, M::step_e());
+        let device = &handles.device;
+        let queue = &handles.queue;
+        let step_bind_group_layout = &handles.step_bind_group_layout;
+        let step_e_pipeline = &handles.step_e_pipeline;
+        let step_h_pipeline = &handles.step_h_pipeline;
+        let sim_meta_buffer = device.create_buffer_init(&wgpu::util::BufferInitDescriptor {
+            label: Some("gpu-side simulation metadata"),
+            contents: to_bytes(&[meta][..]),
+            usage: wgpu::BufferUsages::STORAGE,
+        });
+
+        let stim_e_buffer = device.create_buffer_init(&wgpu::util::BufferInitDescriptor {
+            label: Some("gpu-side stimulus e field"),
+            contents: to_bytes(stim_cpu_e),
+            usage: wgpu::BufferUsages::STORAGE
+        });
+
+        let stim_h_buffer = device.create_buffer_init(&wgpu::util::BufferInitDescriptor {
+            label: Some("gpu-side stimulus h field"),
+            contents: to_bytes(stim_cpu_h),
+            usage: wgpu::BufferUsages::STORAGE
+        });
+
+        let mat_buffer = device.create_buffer_init(&wgpu::util::BufferInitDescriptor {
+            label: Some("gpu-side materials matrix"),
+            contents: to_bytes(mat),
+            // Can be used by the GPU and copied back to the CPU
+            usage: wgpu::BufferUsages::STORAGE,
+        });
+
+        let e_field_buffer = device.create_buffer_init(&wgpu::util::BufferInitDescriptor {
+            label: Some("gpu-side in/out e field"),
+            contents: to_bytes(e),
+            usage: wgpu::BufferUsages::STORAGE.union(wgpu::BufferUsages::COPY_SRC),
+        });
+
+        let h_field_buffer = device.create_buffer_init(&wgpu::util::BufferInitDescriptor {
+            label: Some("gpu-side in/out h field"),
+            contents: to_bytes(h),
+            // Can be used by the GPU and copied back to the CPU
+            usage: wgpu::BufferUsages::STORAGE.union(wgpu::BufferUsages::COPY_SRC),
+        });
+
+        let m_field_buffer = device.create_buffer_init(&wgpu::util::BufferInitDescriptor {
+            label: Some("gpu-side in/out m field"),
+            contents: to_bytes(m),
+            // Can be used by the GPU and copied back to the CPU
+            usage: wgpu::BufferUsages::STORAGE.union(wgpu::BufferUsages::COPY_SRC),
+        });
+
+        let e_readback_buffer = device.create_buffer(&wgpu::BufferDescriptor {
+            label: Some("cpu-side copy of e output buffer"),
+            size: field_bytes as wgpu::BufferAddress,
+            // Can be read to the CPU, and can be copied from the shader's storage buffer
+            usage: wgpu::BufferUsages::MAP_READ | wgpu::BufferUsages::COPY_DST,
+            mapped_at_creation: false,
+        });
+
+        let h_readback_buffer = device.create_buffer(&wgpu::BufferDescriptor {
+            label: Some("cpu-side copy of h output buffer"),
+            size: field_bytes as wgpu::BufferAddress,
+            // Can be read to the CPU, and can be copied from the shader's storage buffer
+            usage: wgpu::BufferUsages::MAP_READ | wgpu::BufferUsages::COPY_DST,
+            mapped_at_creation: false,
+        });
+
+        let m_readback_buffer = device.create_buffer(&wgpu::BufferDescriptor {
+            label: Some("cpu-side copy of m output buffer"),
+            size: field_bytes as wgpu::BufferAddress,
+            // Can be read to the CPU, and can be copied from the shader's storage buffer
+            usage: wgpu::BufferUsages::MAP_READ | wgpu::BufferUsages::COPY_DST,
+            mapped_at_creation: false,
+        });
+
+        let bind_group = device.create_bind_group(&wgpu::BindGroupDescriptor {
+            label: None,
+            layout: &step_bind_group_layout,
+            entries: &[
+                wgpu::BindGroupEntry {
+                    binding: 0,
+                    resource: sim_meta_buffer.as_entire_binding(),
+                },
+                wgpu::BindGroupEntry {
+                    binding: 1,
+                    resource: stim_e_buffer.as_entire_binding(),
+                },
+                wgpu::BindGroupEntry {
+                    binding: 2,
+                    resource: stim_h_buffer.as_entire_binding(),
+                },
+                wgpu::BindGroupEntry {
+                    binding: 3,
+                    resource: mat_buffer.as_entire_binding(),
+                },
+                wgpu::BindGroupEntry {
+                    binding: 4,
+                    resource: e_field_buffer.as_entire_binding(),
+                },
+                wgpu::BindGroupEntry {
+                    binding: 5,
+                    resource: h_field_buffer.as_entire_binding(),
+                },
+                wgpu::BindGroupEntry {
+                    binding: 6,
+                    resource: m_field_buffer.as_entire_binding(),
+                },
+            ],
+        });
+
+        let workgroups = ((meta.dim.x()+3) / 4, (meta.dim.y()+3) / 4, (meta.dim.z()+3) / 4);
+
+        let mut encoder =
+            device.create_command_encoder(&wgpu::CommandEncoderDescriptor { label: None });
+
+        for _ in 0..num_steps {
+            {
+                let mut cpass = encoder.begin_compute_pass(&wgpu::ComputePassDescriptor { label: None });
+                cpass.set_bind_group(0, &bind_group, &[]);
+                cpass.set_pipeline(&step_e_pipeline);
+                cpass.dispatch(workgroups.0, workgroups.1, workgroups.2);
+            }
+
+            {
+                let mut cpass = encoder.begin_compute_pass(&wgpu::ComputePassDescriptor { label: None });
+                cpass.set_bind_group(0, &bind_group, &[]);
+                cpass.set_pipeline(&step_h_pipeline);
+                cpass.dispatch(workgroups.0, workgroups.1, workgroups.2);
+            }
+        }
+
+        encoder.copy_buffer_to_buffer(
+            &e_field_buffer,
+            0,
+            &e_readback_buffer,
+            0,
+            field_bytes as u64,
+        );
+
+        encoder.copy_buffer_to_buffer(
+            &h_field_buffer,
+            0,
+            &h_readback_buffer,
+            0,
+            field_bytes as u64,
+        );
+
+        encoder.copy_buffer_to_buffer(
+            &m_field_buffer,
+            0,
+            &m_readback_buffer,
+            0,
+            field_bytes as u64,
+        );
+
+        queue.submit(Some(encoder.finish()));
+
+        let e_readback_slice = e_readback_buffer.slice(..);
+        let e_readback_future = e_readback_slice.map_async(wgpu::MapMode::Read).then(|_| async {
+                e.copy_from_slice(unsafe {
+                    from_bytes(e_readback_slice.get_mapped_range().as_ref())
+                });
+                e_readback_buffer.unmap();
+            });
+
+        let h_readback_slice = h_readback_buffer.slice(..);
+        let h_readback_future = h_readback_slice.map_async(wgpu::MapMode::Read).then(|_| async {
+            h.copy_from_slice(unsafe {
+                from_bytes(h_readback_slice.get_mapped_range().as_ref())
+            });
+            h_readback_buffer.unmap();
+        });
+
+        let m_readback_slice = m_readback_buffer.slice(..);
+        let m_readback_future = m_readback_slice.map_async(wgpu::MapMode::Read).then(|_| async {
+            m.copy_from_slice(unsafe {
+                from_bytes(m_readback_slice.get_mapped_range().as_ref())
+            });
+            m_readback_buffer.unmap();
+        });
+
+        device.poll(wgpu::Maintain::Wait);
+
+        futures::executor::block_on(futures::future::join(
+            e_readback_future, futures::future::join(
+                h_readback_future, m_readback_future)));
+    }
+}
+
+/// Convert an arbitrary slice into a byte slice
+fn to_bytes<T>(slice: &[T]) -> &[u8] {
+    unsafe {
+        std::slice::from_raw_parts(slice.as_ptr() as *const u8, slice.len() * std::mem::size_of::<T>())
+    }
+}
+
+/// Convert a byte slice into a T slice
+unsafe fn from_bytes<T>(slice: &[u8]) -> &[T] {
+    let elem_size = std::mem::size_of::<T>();
+    let new_len = slice.len() / elem_size;
+    assert_eq!(new_len * elem_size, slice.len());
+    std::slice::from_raw_parts(slice.as_ptr() as *const T, new_len)
+}
+
+/// Loads the shader
+fn get_shader() -> wgpu::ShaderModuleDescriptorSpirV<'static> {
+    let data = spirv_backend_runner::spirv_module();
+    let spirv = Cow::Owned(wgpu::util::make_spirv_raw(&data).into_owned());
+
+    wgpu::ShaderModuleDescriptorSpirV {
+        label: None,
+        source: spirv,
+    }
+}
+
+async fn open_device(max_buf_size: u64) -> (wgpu::Device, wgpu::Queue) {
+    // based on rust-gpu/examples/runners/wgpu/src/compute.rs:start_internal
+    let instance = wgpu::Instance::new(wgpu::Backends::PRIMARY);
+    info!("open_device: got instance");
+    let adapter = instance
+        .request_adapter(&wgpu::RequestAdapterOptions {
+            power_preference: wgpu::PowerPreference::HighPerformance,
+            force_fallback_adapter: false,
+            compatible_surface: None,
+        })
+        .await
+        .expect("Failed to find an appropriate adapter");
+    info!("open_device: got adapter");
+
+    // XXX not all adapters will support non-default limits, and it could
+    // cause perf degradations even on the ones that do. May want to consider
+    // folding some buffers together to avoid this.
+    let mut limits = wgpu::Limits::default();
+    //limits.max_bind_groups = 5;
+    //limits.max_dynamic_storage_buffers_per_pipeline_layout = 5;
+    limits.max_storage_buffers_per_shader_stage = 7;
+    //limits.max_storage_buffer_binding_size = 128 MiB.
+    //limits.max_storage_buffer_binding_size = 1024 * (1 << 20);
+    limits.max_storage_buffer_binding_size = max_buf_size.try_into().unwrap();
+
+    let (device, queue) = adapter
+        .request_device(
+            &wgpu::DeviceDescriptor {
+                label: None,
+                features: wgpu::Features::SPIRV_SHADER_PASSTHROUGH,
+                limits,
+            },
+            None,
+        )
+        .await
+        .expect("Failed to create device");
+    info!("open_device: got device");
+    (device, queue)
+}
+
+fn make_pipelines(
+    device: &wgpu::Device,
+    shader_module: &wgpu::ShaderModule,
+    entry_step_h: &'static str,
+    entry_step_e: &'static str
+) -> (
+    wgpu::BindGroupLayout, wgpu::ComputePipeline, wgpu::ComputePipeline
+) {
+    let bind_group_layout = device.create_bind_group_layout(&wgpu::BindGroupLayoutDescriptor {
+        label: None,
+        entries: &[
+            wgpu::BindGroupLayoutEntry {
+                // meta
+                binding: 0,
+                count: None,
+                visibility: wgpu::ShaderStages::COMPUTE,
+                ty: wgpu::BindingType::Buffer {
+                    has_dynamic_offset: false,
+                    min_binding_size: Some(NonZeroU64::new(1).unwrap()),
+                    ty: wgpu::BufferBindingType::Storage { read_only: true },
+                },
+            },
+            wgpu::BindGroupLayoutEntry {
+                // stimulus(e)
+                binding: 1,
+                count: None,
+                visibility: wgpu::ShaderStages::COMPUTE,
+                ty: wgpu::BindingType::Buffer {
+                    has_dynamic_offset: false,
+                    min_binding_size: Some(NonZeroU64::new(1).unwrap()),
+                    ty: wgpu::BufferBindingType::Storage { read_only: true },
+                },
+            },
+            wgpu::BindGroupLayoutEntry {
+                // stimulus(h)
+                binding: 2,
+                count: None,
+                visibility: wgpu::ShaderStages::COMPUTE,
+                ty: wgpu::BindingType::Buffer {
+                    has_dynamic_offset: false,
+                    min_binding_size: Some(NonZeroU64::new(1).unwrap()),
+                    ty: wgpu::BufferBindingType::Storage { read_only: true },
+                },
+            },
+            wgpu::BindGroupLayoutEntry {
+                // materials
+                binding: 3,
+                count: None,
+                visibility: wgpu::ShaderStages::COMPUTE,
+                ty: wgpu::BindingType::Buffer {
+                    has_dynamic_offset: false,
+                    min_binding_size: Some(NonZeroU64::new(1).unwrap()),
+                    ty: wgpu::BufferBindingType::Storage { read_only: true },
+                },
+            },
+            wgpu::BindGroupLayoutEntry {
+                // e field
+                binding: 4,
+                count: None,
+                visibility: wgpu::ShaderStages::COMPUTE,
+                ty: wgpu::BindingType::Buffer {
+                    has_dynamic_offset: false,
+                    min_binding_size: Some(NonZeroU64::new(1).unwrap()),
+                    ty: wgpu::BufferBindingType::Storage { read_only: false },
+                },
+            },
+            wgpu::BindGroupLayoutEntry {
+                // h field
+                binding: 5,
+                count: None,
+                visibility: wgpu::ShaderStages::COMPUTE,
+                ty: wgpu::BindingType::Buffer {
+                    has_dynamic_offset: false,
+                    min_binding_size: Some(NonZeroU64::new(1).unwrap()),
+                    ty: wgpu::BufferBindingType::Storage { read_only: false },
+                },
+            },
+            wgpu::BindGroupLayoutEntry {
+                // m field
+                binding: 6,
+                count: None,
+                visibility: wgpu::ShaderStages::COMPUTE,
+                ty: wgpu::BindingType::Buffer {
+                    has_dynamic_offset: false,
+                    min_binding_size: Some(NonZeroU64::new(1).unwrap()),
+                    ty: wgpu::BufferBindingType::Storage { read_only: false },
+                },
+            },
+        ],
+    });
+
+    let pipeline_layout = device.create_pipeline_layout(&wgpu::PipelineLayoutDescriptor {
+        label: None,
+        bind_group_layouts: &[&bind_group_layout],
+        push_constant_ranges: &[],
+    });
+
+    let compute_step_h_pipeline = device.create_compute_pipeline(&wgpu::ComputePipelineDescriptor {
+        label: None,
+        layout: Some(&pipeline_layout),
+        module: shader_module,
+        entry_point: entry_step_h,
+    });
+
+    let compute_step_e_pipeline = device.create_compute_pipeline(&wgpu::ComputePipelineDescriptor {
+        label: None,
+        layout: Some(&pipeline_layout),
+        module: shader_module,
+        entry_point: entry_step_e,
+    });
+
+    (bind_group_layout, compute_step_h_pipeline, compute_step_e_pipeline)
+}
+

crates/coremem/src/stim.rs

@@ -1,11 +1,24 @@
-use crate::real::*;
-use crate::geom::{Meters, Region, Vec3};
+use crate::real::Real as _;
+use crate::cross::vec::Vec3;
+use crate::geom::{HasCrossSection, Meters, Region};
 use rand;
 
-type Fields = (Vec3<f32>, Vec3<f32>);
+/// field densities
+#[derive(Clone, Copy, Debug, Default, PartialEq)]
+pub struct Fields {
+    pub e: Vec3<f32>,
+    pub h: Vec3<f32>,
+}
+
+/// field magnitude densities (really, signed magnitude)
+#[derive(Clone, Copy, Debug, Default, PartialEq)]
+pub struct FieldMags {
+    pub e: f32,
+    pub h: f32,
+}
 
 pub trait AbstractStimulus: Sync {
-    // TODO: might be cleaner to return some `Fields` type instead of a tuple
+    // TODO: using floats for time and position is not great.
     /// Return the (E, H) field which should be added PER-SECOND to the provided position/time.
     fn at(&self, t_sec: f32, pos: Meters) -> Fields;
 }
@@ -18,13 +31,13 @@ pub trait AbstractStimulus: Sync {
 
 impl<T: AbstractStimulus> AbstractStimulus for Vec<T> {
     fn at(&self, t_sec: f32, pos: Meters) -> Fields {
-        let (mut e, mut h) = Fields::default();
+        let mut fields = Fields::default();
         for i in self {
-            let (de, dh) = i.at(t_sec, pos);
-            e += de;
-            h += dh;
+            let Fields { e, h } = i.at(t_sec, pos);
+            fields.e += e;
+            fields.h += h;
         }
-        (e, h)
+        fields
     }
 }
 
@@ -63,7 +76,7 @@ impl AbstractStimulus for UniformStimulus {
 impl TimeVarying for UniformStimulus {}
 impl TimeVarying3 for UniformStimulus {
     fn at(&self, _t_sec: f32) -> Fields {
-        (self.e, self.h)
+        Fields { e: self.e, h: self.h }
     }
 }
 
@@ -99,7 +112,10 @@ impl RngStimulus {
 
 impl AbstractStimulus for RngStimulus {
     fn at(&self, t_sec: f32, pos: Meters) -> Fields {
-        (self.gen(t_sec, pos, self.e_scale, 0), self.gen(t_sec, pos, self.h_scale, 0x7de3))
+        Fields {
+            e: self.gen(t_sec, pos, self.e_scale, 0),
+            h: self.gen(t_sec, pos, self.h_scale, 0x7de3)
+        }
     }
 }
 
@@ -134,25 +150,25 @@ impl<R: Region + Sync, T: TimeVarying3 + Sync> AbstractStimulus for Stimulus<R,
 pub struct CurlStimulus<R, T> {
     region: R,
     stim: T,
-    center: Meters,
-    axis: Meters,
 }
 
 impl<R, T> CurlStimulus<R, T> {
-    pub fn new(region: R, stim: T, center: Meters, axis: Meters) -> Self {
-        Self { region, stim, center, axis }
+    pub fn new(region: R, stim: T) -> Self {
+        Self { region, stim }
     }
 }
 
-impl<R: Region + Sync, T: TimeVarying1 + Sync> AbstractStimulus for CurlStimulus<R, T> {
+impl<R: Region + HasCrossSection, T: TimeVarying1 + Sync> AbstractStimulus for CurlStimulus<R, T> {
     fn at(&self, t_sec: f32, pos: Meters) -> Fields {
         if self.region.contains(pos) {
-            let (amt_e, amt_h) = self.stim.at(t_sec);
-            let from_center_to_point = *pos - *self.center;
-            let rotational = from_center_to_point.cross(*self.axis);
-            let impulse_e = rotational.with_mag(amt_e.cast());
-            let impulse_h = rotational.with_mag(amt_h.cast());
-            (impulse_e, impulse_h)
+            let FieldMags { e, h } = self.stim.at(t_sec);
+            let rotational = self.region.cross_section_normal(pos).norm();
+            let impulse_e = rotational.with_mag(e.cast()).unwrap_or_default();
+            let impulse_h = rotational.with_mag(h.cast()).unwrap_or_default();
+            Fields {
+                e: impulse_e,
+                h: impulse_h,
+            }
         } else {
             Fields::default()
         }
@@ -170,7 +186,7 @@ pub trait TimeVarying: Sized {
 
 pub trait TimeVarying1: TimeVarying {
     /// Retrieve the (E, H) impulse to apply PER-SECOND at the provided time (in seconds).
-    fn at(&self, t_sec: f32) -> (f32, f32);
+    fn at(&self, t_sec: f32) -> FieldMags;
 }
 
 pub trait TimeVarying3: TimeVarying {
@@ -181,8 +197,8 @@ pub trait TimeVarying3: TimeVarying {
 // assumed to represent the E field
 impl TimeVarying for f32 {}
 impl TimeVarying1 for f32 {
-    fn at(&self, _t_sec: f32) -> (f32, f32) {
-        (*self, 0.0)
+    fn at(&self, _t_sec: f32) -> FieldMags {
+        FieldMags { e: *self, h: 0.0 }
     }
 }
 
@@ -225,20 +241,20 @@ impl<A> Sinusoid<A> {
 impl<A> TimeVarying for Sinusoid<A> {}
 
 impl TimeVarying1 for Sinusoid1 {
-    fn at(&self, t_sec: f32) -> (f32, f32) {
-        (
-            self.amp * (t_sec * self.omega).sin(),
-            0.0,
-        )
+    fn at(&self, t_sec: f32) -> FieldMags {
+        FieldMags {
+            e: self.amp * (t_sec * self.omega).sin(),
+            h: 0.0,
+        }
     }
 }
 
 impl TimeVarying3 for Sinusoid3 {
     fn at(&self, t_sec: f32) -> Fields {
-        (
-            self.amp * (t_sec * self.omega).sin(),
-            Vec3::zero(),
-        )
+        Fields {
+            e: self.amp * (t_sec * self.omega).sin(),
+            h: Vec3::zero(),
+        }
     }
 }
 
@@ -267,20 +283,20 @@ impl<A> Exp<A> {
 impl<A> TimeVarying for Exp<A> {}
 
 impl TimeVarying1 for Exp1 {
-    fn at(&self, t_sec: f32) -> (f32, f32) {
-        (
-            self.amp * (t_sec * -self.tau).exp(),
-            0.0,
-        )
+    fn at(&self, t_sec: f32) -> FieldMags {
+        FieldMags {
+            e: self.amp * (t_sec * -self.tau).exp(),
+            h: 0.0,
+        }
     }
 }
 
 impl TimeVarying3 for Exp3 {
     fn at(&self, t_sec: f32) -> Fields {
-        (
-            self.amp * (t_sec * -self.tau).exp(),
-            Vec3::zero(),
-        )
+        Fields {
+            e: self.amp * (t_sec * -self.tau).exp(),
+            h: Vec3::zero(),
+        }
     }
 }
 
@@ -299,7 +315,7 @@ impl<T> Gated<T> {
 
 impl<T> TimeVarying for Gated<T> {}
 impl<T: TimeVarying1> TimeVarying1 for Gated<T> {
-    fn at(&self, t_sec: f32) -> (f32, f32) {
+    fn at(&self, t_sec: f32) -> FieldMags {
         if (self.start..self.end).contains(&t_sec) {
             self.inner.at(t_sec)
         } else {
@@ -332,7 +348,7 @@ impl<T> Shifted<T> {
 
 impl<T> TimeVarying for Shifted<T> {}
 impl<T: TimeVarying1> TimeVarying1 for Shifted<T> {
-    fn at(&self, t_sec: f32) -> (f32, f32) {
+    fn at(&self, t_sec: f32) -> FieldMags {
         self.inner.at(t_sec - self.start_at)
     }
 }
@@ -346,15 +362,16 @@ impl<T: TimeVarying3> TimeVarying3 for Shifted<T> {
 #[cfg(test)]
 mod test {
     use super::*;
+    use serde::{Serialize, Deserialize};
 
     macro_rules! assert_approx_eq {
         ($x:expr, $e:expr, $h:expr) => {
             let x = $x;
             let e = $e;
             let h = $h;
-            let diff_e = (x.0 - e).mag();
+            let diff_e = (x.e - e).mag();
             assert!(diff_e <= 0.001, "{:?} != {:?}", x, e);
-            let diff_h = (x.1 - h).mag();
+            let diff_h = (x.h - h).mag();
             assert!(diff_h <= 0.001, "{:?} != {:?}", x, h);
         }
     }
@@ -362,7 +379,7 @@ mod test {
     #[test]
     fn sinusoid3() {
         let s = Sinusoid3::new(Vec3::new(10.0, 1.0, -100.0), 1000.0);
-        assert_eq!(s.at(0.0), (Vec3::zero(), Vec3::zero()));
+        assert_eq!(s.at(0.0), Fields::default());
         assert_approx_eq!(s.at(0.00025),
             Vec3::new(10.0, 1.0, -100.0), Vec3::zero()
         );
@@ -373,9 +390,67 @@ mod test {
     #[test]
     fn sinusoid3_from_wavelength() {
         let s = Sinusoid3::from_wavelength(Vec3::new(10.0, 1.0, -100.0), 0.001);
-        assert_eq!(s.at(0.0), (Vec3::zero(), Vec3::zero()));
+        assert_eq!(s.at(0.0), Fields::default());
         assert_approx_eq!(s.at(0.00025), Vec3::new(10.0, 1.0, -100.0), Vec3::zero());
         assert_approx_eq!(s.at(0.00050), Vec3::zero(), Vec3::zero());
         assert_approx_eq!(s.at(0.00075), Vec3::new(-10.0, -1.0, 100.0), Vec3::zero());
     }
+
+    #[derive(Clone, Serialize, Deserialize)]
+    struct MockRegion {
+        normal: Vec3<f32>,
+    }
+    impl HasCrossSection for MockRegion {
+        fn cross_section_normal(&self, _p: Meters) -> Vec3<f32> {
+            self.normal
+        }
+    }
+    #[typetag::serde]
+    impl Region for MockRegion {
+        fn contains(&self, _p: Meters) -> bool {
+            true
+        }
+    }
+
+    #[test]
+    fn curl_stimulus_trivial() {
+        let region = MockRegion {
+            normal: Vec3::new(1.0, 0.0, 0.0)
+        };
+        // stim acts as e: 1.0, h: 0.0
+        let stim = CurlStimulus::new(region, 1.0);
+        assert_eq!(stim.at(0.0, Meters::new(0.0, 0.0, 0.0)), Fields {
+            e: Vec3::new(1.0, 0.0, 0.0),
+            h: Vec3::new(0.0, 0.0, 0.0),
+        });
+    }
+
+    #[test]
+    fn curl_stimulus_scales_right() {
+        let region = MockRegion {
+            // the magnitude of this normal shouldn't influence the curl.
+            // though, maybe it would be more useful if it *did*?
+            normal: Vec3::new(5.0, 0.0, 0.0)
+        };
+        // stim acts as e: -3.0, h: 0.0
+        let stim = CurlStimulus::new(region, -3.0);
+        assert_eq!(stim.at(0.0, Meters::new(0.0, 0.0, 0.0)), Fields {
+            e: Vec3::new(-3.0, 0.0, 0.0),
+            h: Vec3::new(0.0, 0.0, 0.0),
+        });
+    }
+
+    #[test]
+    fn curl_stimulus_multi_axis() {
+        let region = MockRegion {
+            normal: Vec3::new(0.0, -1.0, 1.0)
+        };
+        // stim acts as e: 2.0, h: 0.0
+        let stim = CurlStimulus::new(region, 2.0);
+        let Fields { e, h: _ } = stim.at(0.0, Meters::new(0.0, 0.0, 0.0));
+        assert!(e.distance(
+                Vec3::new(0.0, -1.0, 1.0).with_mag(2.0).unwrap()
+            ) < 1e-6
+        );
+    }
 }

crates/coremem/src/util/mod.rs

@@ -1 +0,0 @@
-pub mod cache;

crates/cross/Cargo.toml

@@ -1,10 +1,11 @@
 [package]
-name = "coremem_types"
-version = "0.1.0"
+name = "coremem_cross"
+version = "0.2.0"
 authors = ["Colin <colin@uninsane.org>"]
 edition = "2021"
 
 [features]
+# some functionality does not compile for the spirv target, so we feature gate these.
 serde = [ "dep:serde" ]
 fmt = []
 

crates/cross/src/compound/enumerated.rs

@@ -66,7 +66,7 @@ where
 }
 
 impl DiscrDispatch<P0> for Discr<P0> {
-    fn dispatch<H: DiscrHandler<P0, O>, O>(&self, h: H) -> O {
+    fn dispatch<H: DiscrHandler<P0, O>, O>(&self, _h: H) -> O {
         unreachable!()
     }
 }
@@ -148,9 +148,12 @@ where
 #[derive(Copy, Clone, Default, PartialEq)]
 pub struct Enum<D, L>(D, L);
 
+/// Users should prefer this type rather than rely on the internal Enum struct implementation.
+#[allow(dead_code)]
 pub type InternallyDiscriminated<Args> = Enum<(), List<Args>>;
 
 impl<P: Peano, L> Enum<(Discr<P>,), L> {
+    #![allow(dead_code)]
     pub fn new<Variants>(v: Variants) -> Self
     where
         Variants: IntoList<List=L>,
@@ -160,6 +163,7 @@ impl<P: Peano, L> Enum<(Discr<P>,), L> {
     }
 }
 impl<L> Enum<(), L> {
+    #![allow(dead_code)]
     pub fn internally_discriminated<Variants>(v: Variants) -> Self
     where
         Variants: IntoList<List=L>,
@@ -168,13 +172,6 @@ impl<L> Enum<(), L> {
     }
 }
 
-impl<D, L> Enum<D, L> {
-    /// use with care. long term, this probably shouldn't stay public.
-    pub fn from_components(discr: D, list: L) -> Self {
-        Self(discr, list)
-    }
-}
-
 pub trait EnumRequirements {
     type NumVariants: Peano;
     fn decode_discr(&self) -> Discr<Self::NumVariants>;
@@ -225,6 +222,7 @@ where
     }
 
     /// invoke the closure on the active variant, passing the variant by mutable reference
+    #[allow(dead_code)]
     pub fn dispatch_mut<'a, F, R>(&'a mut self, f: F) -> R
     where
         DispatchIndexable<&'a mut L, F>: DiscrHandler<<Self as EnumRequirements>::NumVariants, R>,

crates/cross/src/compound/list/flat.rs

@@ -127,15 +127,20 @@ pub trait IntoList {
     fn into_list(self) -> Self::List;
 }
 
-pub type List<Args> = <Args as IntoList>::List;
-pub type List1<E0> = Node<E0, Null>;
-pub type List2<E0, E1> = Node<E0, List1<E1>>;
-pub type List3<E0, E1, E2> = Node<E0, List2<E1, E2>>;
-pub type List4<E0, E1, E2, E3> = Node<E0, List3<E1, E2, E3>>;
-
 impl<H> Meta for Node<H, Null> {
     type Length = P1;
 }
 impl<H, T: Meta> Meta for Node<H, T> {
     type Length = PNext<T::Length>;
 }
+
+/// these are exported for the convenience of potential consumers: not needed internally
+pub(crate) mod exports {
+    #![allow(dead_code)]
+    use super::{IntoList, Node, Null};
+    pub type List<Args> = <Args as IntoList>::List;
+    pub type List1<E0> = Node<E0, Null>;
+    pub type List2<E0, E1> = Node<E0, List1<E1>>;
+    pub type List3<E0, E1, E2> = Node<E0, List2<E1, E2>>;
+    pub type List4<E0, E1, E2, E3> = Node<E0, List3<E1, E2, E3>>;
+}

crates/cross/src/compound/list/mod.rs

@@ -1,11 +1,12 @@
 use crate::compound::peano::{Peano, PeanoNonZero};
 
 mod flat;
-mod linked;
-mod tuple_consumer;
+// mod linked;
+// mod tuple_consumer;
 // pub use tuple_consumer::*;
 // pub use linked::*;
-pub use flat::*;
+pub use flat::IntoList;
+pub use flat::exports::*;
 
 pub trait Indexable<P: Peano> {
     type Element;

crates/cross/src/compound/mod.rs

@@ -1,3 +1,6 @@
 pub mod enumerated;
 pub mod list;
+mod optional;
 pub mod peano;
+
+pub use optional::Optional;

crates/cross/src/compound/optional.rs

@@ -0,0 +1,86 @@
+#[cfg(feature = "serde")]
+use serde::{Serialize, Deserialize};
+
+/// This is a spirv-compatible option type.
+/// The native rust Option type produces invalid spirv due to its enum nature; this custom option
+/// type creates code which will actually compile.
+#[cfg_attr(feature = "serde", derive(Deserialize, Serialize))]
+#[cfg_attr(feature = "fmt", derive(Debug))]
+#[derive(Copy, Clone, PartialEq)]
+pub struct Optional<T> {
+    // XXX: not a bool, because: "entrypoint parameter cannot contain a boolean"
+    present: u8,
+    data: T,
+}
+
+impl<T> Optional<T> {
+    pub fn some(data: T) -> Self {
+        Self {
+            present: 1,
+            data,
+        }
+    }
+
+    pub fn explicit_none(data: T) -> Self {
+        Self {
+            present: 0,
+            data,
+        }
+    }
+
+    pub fn is_some(self) -> bool {
+        self.present != 0
+    }
+
+    pub fn unwrap(self) -> T {
+        debug_assert!(self.present != 0);
+        self.data
+    }
+
+    pub fn map<U: Default, F: FnOnce(T) -> U>(self, f: F) -> Optional<U> {
+        self.and_then(|inner| Optional::some(f(inner)))
+    }
+
+    pub fn and_then<U: Default, F: FnOnce(T) -> Optional<U>>(self, f: F) -> Optional<U> {
+        if self.present != 0 {
+            f(self.data)
+        } else {
+            Optional::none()
+        }
+    }
+
+    pub fn unwrap_or(self, default: T) -> T {
+        if self.present != 0 {
+            self.data
+        } else {
+            default
+        }
+    }
+}
+
+impl<T: Default> Optional<T> {
+    pub fn none() -> Self {
+        Self::explicit_none(Default::default())
+    }
+
+    pub fn unwrap_or_default(self) -> T {
+        self.unwrap_or(Default::default())
+    }
+}
+
+impl<T: Default> Default for Optional<T> {
+    fn default() -> Self {
+        Self::none()
+    }
+}
+
+impl<T0: Default, T1: Default> Optional<(T0, T1)> {
+    pub fn flatten((f0, f1): (Optional<T0>, Optional<T1>)) -> Self {
+        if f0.present != 0 && f1.present != 0 {
+            Optional::some((f0.data, f1.data))
+        } else {
+            Optional::none()
+        }
+    }
+}
+

crates/cross/src/compound/peano.rs

@@ -0,0 +1,68 @@
+//! Peano numbers (also known as Church numerals) are type-level natural numbers.
+//! each non-zero Peano number is defined as the unique successor of a previous Peano number.
+//!
+//! - given some Peano number I, we can derive its successor by S=PNext<I>.
+//! - given a Peano number PNext<I>, we can define its predecessor as I=P.
+//! - the base Peano number, which represents 0, is `P0`.
+//! - the `Peano` trait exposes alternative syntaxes for these: `P::Next` and `P::Prev`,
+//!   however the type system can reason less about these (they're just a convenience).
+//!
+//! the primary use of Peano numbers is to allow types to specialize on a specific natural number
+//! out of some larger set of natural numbers. e.g. one might have a `struct List<Length: Peano>`
+//! to allow constructing a list of compile-time constant length.
+//!
+//! this whole module will hopefully be obsoleted as Rust's type-level integers become more
+//! capable, but in 2022 Peano numbers enable more operations (arithmetic, specialization) than type-level integers.
+
+#[derive(Copy, Clone, Default, PartialEq)]
+pub struct PNext<P>(P);
+#[derive(Copy, Clone, Default, PartialEq)]
+pub struct P0;
+pub type P1 = PNext<P0>;
+
+/// these are exported for the convenience of potential consumers: not needed internally
+mod exports {
+    #![allow(dead_code)]
+    use super::{P1, PNext};
+    pub type P2 = PNext<P1>;
+    pub type P3 = PNext<P2>;
+    pub type P4 = PNext<P3>;
+    pub type P5 = PNext<P4>;
+    pub type P6 = PNext<P5>;
+    pub type P7 = PNext<P6>;
+}
+pub use exports::*;
+
+pub trait Peano: Copy + Clone + Default + PartialEq {
+    type Next: PeanoNonZero;
+    type PrevOrZero: Peano;
+    /// always set to ()
+    /// this exists to allow Peano numbers to be used as struct parameters without PhantomData
+    type Unit: Copy + Clone + Default + PartialEq;
+    const VALUE: u32;
+}
+
+pub trait PeanoNonZero: Peano {
+    type Prev: Peano;
+}
+
+impl Peano for P0 {
+    type Next = P1;
+    type PrevOrZero = P0;
+    type Unit = ();
+    const VALUE: u32 = 0;
+}
+impl<P: Peano> Peano for PNext<P> {
+    type Next = PNext<PNext<P>>;
+    type PrevOrZero = P;
+    type Unit = ();
+    const VALUE: u32 = 1 + P::VALUE;
+}
+impl<P: Peano> PeanoNonZero for PNext<P> {
+    type Prev = P;
+}
+
+// A: LessThan<B> is satisfied only if A is strictly less than B.
+pub trait LessThan<P: Peano> { }
+impl<P: Peano> LessThan<PNext<P>> for P { }
+

crates/cross/src/dim.rs

@@ -0,0 +1,66 @@
+use core::ops::{Index, IndexMut};
+
+use crate::vec::Vec3u;
+
+/// use this to wrap a flat region of memory into something which can be indexed by coordinates in
+/// 3d space.
+pub struct DimensionedSlice<T> {
+    dim: Vec3u,
+    items: T,
+}
+
+impl<T> DimensionedSlice<T> {
+    pub fn new(dim: Vec3u, items: T) -> Self {
+        Self { dim, items }
+    }
+}
+
+impl<'a, T: Index<usize> + ?Sized> Index<Vec3u> for DimensionedSlice<&'a T> {
+    type Output=T::Output;
+
+    fn index(&self, idx: Vec3u) -> &Self::Output {
+        let idx = index(idx, self.dim);
+        &self.items[idx]
+    }
+}
+
+impl<'a, T: Index<usize> + ?Sized> Index<Vec3u> for DimensionedSlice<&'a mut T> {
+    type Output=T::Output;
+
+    fn index(&self, idx: Vec3u) -> &Self::Output {
+        let idx = index(idx, self.dim);
+        &self.items[idx]
+    }
+}
+
+impl<'a, T: IndexMut<usize> + ?Sized> IndexMut<Vec3u> for DimensionedSlice<&'a mut T> {
+    fn index_mut(&mut self, idx: Vec3u) -> &mut Self::Output {
+        let idx = index(idx, self.dim);
+        &mut self.items[idx]
+    }
+}
+
+fn index(loc: Vec3u, dim: Vec3u) -> usize {
+    ((loc.z()*dim.y() + loc.y())*dim.x() + loc.x()) as usize
+}
+
+
+#[cfg(test)]
+mod test {
+    use super::*;
+    #[test]
+    fn test_index() {
+        let dim = Vec3u::new(2, 3, 7);
+        assert_eq!(index(Vec3u::new(0, 0, 0), dim), 0);
+        assert_eq!(index(Vec3u::new(1, 0, 0), dim), 1);
+        assert_eq!(index(Vec3u::new(0, 1, 0), dim), 2);
+        assert_eq!(index(Vec3u::new(1, 1, 0), dim), 3);
+        assert_eq!(index(Vec3u::new(0, 2, 0), dim), 4);
+        assert_eq!(index(Vec3u::new(0, 0, 1), dim), 6);
+        assert_eq!(index(Vec3u::new(1, 0, 1), dim), 7);
+        assert_eq!(index(Vec3u::new(0, 1, 1), dim), 8);
+        assert_eq!(index(Vec3u::new(1, 2, 1), dim), 11);
+        assert_eq!(index(Vec3u::new(1, 2, 2), dim), 17);
+    }
+}
+

crates/cross/src/lib.rs

@@ -0,0 +1,11 @@
+#![no_std]
+#![feature(core_intrinsics)]
+
+pub mod compound;
+pub mod dim;
+pub mod mat;
+pub mod real;
+pub mod step;
+pub mod vec;
+// private because `vec` re-exports to important vecu constructs
+mod vecu;

crates/cross/src/mat/compound.rs

@@ -26,12 +26,6 @@ pub struct DiscrMat<Mats>(Enum<(), Mats>);
 pub type DiscrMat2<M0, M1> = DiscrMat<List2<M0, M1>>;
 pub type DiscrMat3<M0, M1, M2> = DiscrMat<List3<M0, M1, M2>>;
 
-impl<Mats> DiscrMat<Mats> {
-    fn new_from_fields(mats: Mats) -> Self {
-        Self(Enum::from_components((), mats))
-    }
-}
-
 impl<Mats: Default> DiscrMat<Mats> {
     fn new<P: Peano>(m: Mats::Element) -> Self
     where
@@ -63,6 +57,24 @@ impl<R: Real, P: Peano, T: Material<R>> VariantHandler<P, T, Vec3<R>> for MoveBV
     }
 }
 
+/// invokes Into<T>::into on any variant
+struct IntoDispatcher;
+impl<P: Peano, T: Into<I>, I> VariantHandler<P, T, I> for IntoDispatcher {
+    fn call(self, v: T) -> I {
+        v.into()
+    }
+}
+
+impl<R, M0, M1> Into<FullyGenericMaterial<R>> for DiscrMat2<M0, M1>
+where
+    M0: DiscriminantCodable<P2> + Into<FullyGenericMaterial<R>> + Copy,
+    M1: Into<FullyGenericMaterial<R>> + Copy,
+{
+    fn into(self) -> FullyGenericMaterial<R> {
+        self.0.dispatch(IntoDispatcher)
+    }
+}
+
 impl<R: Real, M0, M1> Material<R> for DiscrMat2<M0, M1>
 where
     M0: DiscriminantCodable<P2> + Material<R> + Copy,
@@ -76,6 +88,17 @@ where
     }
 }
 
+impl<R, M0, M1, M2> Into<FullyGenericMaterial<R>> for DiscrMat3<M0, M1, M2>
+where
+    M0: DiscriminantCodable<P3> + Into<FullyGenericMaterial<R>> + Copy,
+    M1: Into<FullyGenericMaterial<R>> + Copy,
+    M2: Into<FullyGenericMaterial<R>> + Copy,
+{
+    fn into(self) -> FullyGenericMaterial<R> {
+        self.0.dispatch(IntoDispatcher)
+    }
+}
+
 impl<R: Real, M0, M1, M2> Material<R> for DiscrMat3<M0, M1, M2>
 where
     M0: DiscriminantCodable<P3> + Material<R> + Copy,
@@ -139,7 +162,7 @@ pub struct GenericMagnetic<R>(DiscrMat3<MBPgram<R>, MHPgram<R>, Vacuum>);
 impl<R: Real> Default for GenericMagnetic<R> {
     fn default() -> Self {
         // N.B.: the default is not the first variant.
-        // we order the variants specifically so that the first one can store the descriminant, but
+        // we order the variants specifically so that the first one can store the discriminant, but
         // we NEED Vacuum to be the default.
         Vacuum.into()
     }
@@ -192,6 +215,11 @@ impl<R: Real> From<MHPgram<R>> for FullyGenericMaterial<R> {
         Self::new(Default::default(), mat.into())
     }
 }
+impl<R: Real> From<Vacuum> for FullyGenericMaterial<R> {
+    fn from(mat: Vacuum) -> Self {
+        Self::new(Default::default(), mat.into())
+    }
+}
 impl<R: Real> From<IsomorphicConductor<R>> for FullyGenericMaterial<R> {
     fn from(mat: IsomorphicConductor<R>) -> Self {
         let mat: AnisomorphicConductor<R> = mat.into();

crates/cross/src/real.rs

@@ -20,11 +20,19 @@ pub trait ToFloat {
     }
 }
 
-#[cfg(feature = "fmt")]
-pub trait RealFeatures: fmt::LowerExp + fmt::Display + fmt::Debug {}
-#[cfg(not(feature = "fmt"))]
+
+#[cfg(all(not(feature = "fmt"), not(feature = "serde")))]
 pub trait RealFeatures {}
 
+#[cfg(all(not(feature = "fmt"), feature = "serde"))]
+pub trait RealFeatures: Serialize + for<'a> Deserialize<'a> {}
+
+#[cfg(all(feature = "fmt", not(feature = "serde")))]
+pub trait RealFeatures: fmt::LowerExp + fmt::Display + fmt::Debug {}
+
+#[cfg(all(feature = "fmt", feature = "serde"))]
+pub trait RealFeatures: fmt::LowerExp + fmt::Display + fmt::Debug + Serialize + for<'a> Deserialize<'a> {}
+
 /// This exists to allow configuration over # of bits (f32 v.s. f64) as well as
 /// constraints.
 pub trait Real:
@@ -99,6 +107,10 @@ pub trait Real:
         self == Self::zero()
     }
 
+    fn inv(self) -> Self {
+        Self::one() / self
+    }
+
     fn zero() -> Self;
     fn one() -> Self;
     fn two() -> Self;
@@ -292,7 +304,7 @@ impl<T: Real> Finite<T> {
     }
     #[cfg(not(feature = "fmt"))]
     fn handle_non_finite(_inner: T) -> ! {
-        panic!("Finite<T> is not finite");
+        panic!(); // expected a finite real
     }
 }
 
@@ -384,7 +396,6 @@ impl<T: ToFloat> ToFloat for Finite<T> {
 impl<T: RealFeatures> RealFeatures for Finite<T> {}
 
 impl<T: Real> Real for Finite<T> {
-    decl_consts!(Self::from_primitive);
     fn from_primitive<P: ToFloat>(p: P) -> Self {
         Self::new(T::from_primitive(p))
     }
@@ -422,6 +433,53 @@ impl<T: Real> Real for Finite<T> {
         let (s, c) = self.0.sin_cos();
         (Self::new(s), Self::new(c))
     }
+
+    // we would ideally use `decl_consts` here, but that produces f64 -> f32 casts for R32 code.
+    fn zero() -> Self {
+        Self::from_primitive(T::zero())
+    }
+    fn one() -> Self {
+        Self::from_primitive(T::one())
+    }
+    fn two() -> Self {
+        Self::from_primitive(T::two())
+    }
+    fn three() -> Self {
+        Self::from_primitive(T::three())
+    }
+    fn ten() -> Self {
+        Self::from_primitive(T::ten())
+    }
+    fn tenth() -> Self {
+        Self::from_primitive(T::tenth())
+    }
+    fn third() -> Self {
+        Self::from_primitive(T::third())
+    }
+    fn half() -> Self {
+        Self::from_primitive(T::half())
+    }
+    fn pi() -> Self {
+        Self::from_primitive(T::pi())
+    }
+    fn two_pi() -> Self {
+        Self::from_primitive(T::two_pi())
+    }
+    fn c() -> Self {
+        Self::from_primitive(T::c())
+    }
+    fn eps0() -> Self {
+        Self::from_primitive(T::eps0())
+    }
+    fn twice_eps0() -> Self {
+        Self::from_primitive(T::twice_eps0())
+    }
+    fn mu0() -> Self {
+        Self::from_primitive(T::mu0())
+    }
+    fn mu0_inv() -> Self {
+        Self::from_primitive(T::mu0_inv())
+    }
 }
 
 impl ToFloat for i32 {

crates/cross/src/step.rs

@@ -0,0 +1,393 @@
+use core::ops::{Index, IndexMut};
+
+use crate::compound::Optional;
+use crate::dim::DimensionedSlice;
+use crate::mat::Material;
+use crate::real::Real;
+use crate::vec::{Vec3, Vec3u};
+
+#[cfg(feature = "serde")]
+use serde::{Serialize, Deserialize};
+
+#[cfg_attr(feature = "serde", derive(Deserialize, Serialize))]
+#[cfg_attr(feature = "fmt", derive(Debug))]
+#[derive(Copy, Clone, Default)]
+pub struct SimMeta<R> {
+    // TODO: make these private?
+    pub dim: Vec3u,
+    pub inv_feature_size: R,
+    pub time_step: R,
+    pub feature_size: R,
+}
+
+impl<R: Real> SimMeta<R> {
+    pub fn new(dim: Vec3u, feature_size: R, time_step: R) -> Self {
+        Self {
+            dim,
+            inv_feature_size: feature_size.inv(),
+            time_step,
+            feature_size
+        }
+    }
+}
+
+impl<R: Copy> SimMeta<R> {
+    pub fn dim(&self) -> Vec3u {
+        self.dim
+    }
+    pub fn inv_feature_size(&self) -> R {
+        self.inv_feature_size
+    }
+    pub fn time_step(&self) -> R {
+        self.time_step
+    }
+    pub fn feature_size(&self) -> R {
+        self.feature_size
+    }
+}
+
+impl<R: Real> SimMeta<R> {
+    pub fn cast<R2: Real>(self) -> SimMeta<R2> {
+        SimMeta {
+            dim: self.dim,
+            inv_feature_size: self.inv_feature_size.cast(),
+            time_step: self.time_step.cast(),
+            feature_size: self.feature_size.cast(),
+        }
+    }
+}
+
+/// Package the field vectors adjacent to some particular location.
+/// Particular those at negative offsets from the midpoint.
+/// This is used in step_e when looking at the H field deltas.
+#[derive(Copy, Clone)]
+pub struct VolumeSampleNeg<R> {
+    pub mid: Vec3<R>,
+    pub xm1: Optional<Vec3<R>>,
+    pub ym1: Optional<Vec3<R>>,
+    pub zm1: Optional<Vec3<R>>,
+}
+
+impl<R: Copy + Default> VolumeSampleNeg<R> {
+    pub fn from_indexable<I: Index<Vec3u, Output=Vec3<R>>>(i: &I, idx: Vec3u) -> Self {
+        VolumeSampleNeg {
+            mid: i[idx],
+            xm1: prev_x(idx).map(|idx| i[idx]),
+            ym1: prev_y(idx).map(|idx| i[idx]),
+            zm1: prev_z(idx).map(|idx| i[idx]),
+        }
+    }
+}
+
+fn prev_x(idx: Vec3u) -> Optional<Vec3u> {
+    match idx.into() {
+        (0, _, _) => Optional::none(),
+        (x, y, z) => Optional::some(Vec3u::new(x-1, y, z)),
+    }
+}
+fn prev_y(idx: Vec3u) -> Optional<Vec3u> {
+    match idx.into() {
+        (_, 0, _) => Optional::none(),
+        (x, y, z) => Optional::some(Vec3u::new(x, y-1, z)),
+    }
+}
+fn prev_z(idx: Vec3u) -> Optional<Vec3u> {
+    match idx.into() {
+        (_, _, 0) => Optional::none(),
+        (x, y, z) => Optional::some(Vec3u::new(x, y, z-1)),
+    }
+}
+
+impl<R: Real> VolumeSampleNeg<R> {
+    /// Calculate the delta in H values amongst this cell and its neighbors (left/up/out)
+    fn delta_h(self) -> FieldDeltas<R> {
+        let mid = self.mid;
+        // let (dfy_dx, dfz_dx) = self.xm1.map(|xm1| {
+        //     (mid.y() - xm1.y(), mid.z() - xm1.z())
+        // }).unwrap_or_default();
+
+        // let (dfx_dy, dfz_dy) = self.ym1.map(|ym1| {
+        //     (mid.x() - ym1.x(), mid.z() - ym1.z())
+        // }).unwrap_or_default();
+
+        // let (dfx_dz, dfy_dz) = self.zm1.map(|zm1| {
+        //     (mid.x() - zm1.x(), mid.y() - zm1.y())
+        // }).unwrap_or_default();
+
+        let (dfy_dx, dfz_dx) = if self.xm1.is_some() {
+            (mid.y() - self.xm1.unwrap().y(), mid.z() - self.xm1.unwrap().z())
+        } else {
+            (R::zero(), R::zero())
+        };
+
+        let (dfx_dy, dfz_dy) = if self.ym1.is_some() {
+            (mid.x() - self.ym1.unwrap().x(), mid.z() - self.ym1.unwrap().z())
+        } else {
+            (R::zero(), R::zero())
+        };
+
+        let (dfx_dz, dfy_dz) = if self.zm1.is_some() {
+            (mid.x() - self.zm1.unwrap().x(), mid.y() - self.zm1.unwrap().y())
+        } else {
+            (R::zero(), R::zero())
+        };
+        
+        FieldDeltas {
+            dfy_dx,
+            dfz_dx,
+            dfx_dy,
+            dfz_dy,
+            dfx_dz,
+            dfy_dz,
+        }
+    }
+}
+
+
+/// Package the field vectors adjacent to some particular location.
+/// Particular those at positive offsets from the midpoint.
+/// This is used in step_h when looking at the E field deltas.
+#[derive(Copy, Clone)]
+pub struct VolumeSamplePos<R> {
+    pub mid: Vec3<R>,
+    pub xp1: Optional<Vec3<R>>,
+    pub yp1: Optional<Vec3<R>>,
+    pub zp1: Optional<Vec3<R>>
+}
+
+impl<R: Copy + Default> VolumeSamplePos<R> {
+    pub fn from_indexable<I: Index<Vec3u, Output=Vec3<R>>>(i: &I, dim: Vec3u, idx: Vec3u) -> Self {
+        VolumeSamplePos {
+            mid: i[idx],
+            xp1: next_x(dim, idx).map(|idx| i[idx]),
+            yp1: next_y(dim, idx).map(|idx| i[idx]),
+            zp1: next_z(dim, idx).map(|idx| i[idx]),
+        }
+    }
+}
+
+fn next_x(dim: Vec3u, idx: Vec3u) -> Optional<Vec3u> {
+    match idx.into() {
+        (x, y, z) if x + 1 < dim.x() => Optional::some(Vec3u::new(x+1, y, z)),
+        _ => Optional::none(),
+    }
+}
+fn next_y(dim: Vec3u, idx: Vec3u) -> Optional<Vec3u> {
+    match idx.into() {
+        (x, y, z) if y + 1 < dim.y() => Optional::some(Vec3u::new(x, y+1, z)),
+        _ => Optional::none(),
+    }
+}
+fn next_z(dim: Vec3u, idx: Vec3u) -> Optional<Vec3u> {
+    match idx.into() {
+        (x, y, z) if z + 1 < dim.z() => Optional::some(Vec3u::new(x, y, z+1)),
+        _ => Optional::none(),
+    }
+}
+
+impl<R: Real> VolumeSamplePos<R> {
+    /// Calculate the delta in E values amongst this cell and its neighbors (right/down/in)
+    fn delta_e(self) -> FieldDeltas<R> {
+        let mid = self.mid;
+        // let (dfy_dx, dfz_dx) = self.xp1.map(|xp1| {
+        //     (xp1.y() - mid.y(), xp1.z() - mid.z())
+        // }).unwrap_or_default();
+
+        // let (dfx_dy, dfz_dy) = self.yp1.map(|yp1| {
+        //     (yp1.x() - mid.x(), yp1.z() - mid.z())
+        // }).unwrap_or_default();
+
+        // let (dfx_dz, dfy_dz) = self.zp1.map(|zp1| {
+        //     (zp1.x() - mid.x(), zp1.y() - mid.y())
+        // }).unwrap_or_default();
+
+        let (dfy_dx, dfz_dx) = if self.xp1.is_some() {
+            (self.xp1.unwrap().y() - mid.y(), self.xp1.unwrap().z() - mid.z())
+        } else {
+            (R::zero(), R::zero())
+        };
+
+        let (dfx_dy, dfz_dy) = if self.yp1.is_some() {
+            (self.yp1.unwrap().x() - mid.x(), self.yp1.unwrap().z() - mid.z())
+        } else {
+            (R::zero(), R::zero())
+        };
+
+        let (dfx_dz, dfy_dz) = if self.zp1.is_some() {
+            (self.zp1.unwrap().x() - mid.x(), self.zp1.unwrap().y() - mid.y())
+        } else {
+            (R::zero(), R::zero())
+        };
+        
+        FieldDeltas {
+            dfy_dx,
+            dfz_dx,
+            dfx_dy,
+            dfz_dy,
+            dfx_dz,
+            dfy_dz,
+        }
+    }
+}
+
+struct FieldDeltas<R> {
+    dfy_dx: R,
+    dfz_dx: R,
+    dfx_dy: R,
+    dfz_dy: R,
+    dfx_dz: R,
+    dfy_dz: R,
+}
+
+impl<R: Real> FieldDeltas<R> {
+    fn nabla(self) -> Vec3<R> {
+        Vec3::new(
+            self.dfz_dy - self.dfy_dz,
+            self.dfx_dz - self.dfz_dx,
+            self.dfy_dx - self.dfx_dy,
+        )
+    }
+}
+
+pub struct StepEContext<'a, R, M> {
+    pub inv_feature_size: R,
+    pub time_step: R,
+    pub stim_e: Vec3<R>,
+    pub mat: &'a M,
+    /// Input field sampled near this location
+    pub in_h: VolumeSampleNeg<R>,
+    pub in_e: Vec3<R>,
+}
+
+impl<'a, R: Real, M: Material<R>> StepEContext<'a, R, M> {
+    pub fn step_flat_view<RM, RF, WF>(
+        meta: SimMeta<R>,
+        mat: &RM,
+        stim_e: &RF,
+        e: &mut WF,
+        h: &RF,
+        idx: Vec3u,
+    )
+    where
+        RM: Index<usize, Output=M> + ?Sized,
+        RF: Index<usize, Output=Vec3<R>> + ?Sized,
+        WF: Index<usize, Output=Vec3<R>> + IndexMut<usize> + ?Sized,
+    {
+        let dim = meta.dim;
+        let stim_e_matrix = DimensionedSlice::new(dim, stim_e);
+        let mat_matrix = DimensionedSlice::new(dim, mat);
+        let mut e_matrix = DimensionedSlice::new(dim, e);
+        let h_matrix = DimensionedSlice::new(dim, h);
+
+        let stim_e = stim_e_matrix[idx];
+        let mat = &mat_matrix[idx];
+        let in_e = e_matrix[idx];
+        let in_h = VolumeSampleNeg::from_indexable(&h_matrix, idx);
+
+        let update_state = StepEContext {
+            inv_feature_size: meta.inv_feature_size,
+            time_step: meta.time_step,
+            stim_e,
+            mat,
+            in_h,
+            in_e,
+        };
+        let new_e = update_state.step_e();
+        e_matrix[idx] = new_e;
+    }
+
+    pub fn step_e(self) -> Vec3<R> {
+        let twice_eps0 = R::twice_eps0();
+        let deltas = self.in_h.delta_h();
+        // \nabla x H
+        let nabla_h = deltas.nabla() * self.inv_feature_size;
+        // $\nabla x H = \epsilon_0 dE/dt + \sigma E$
+        // no-conductivity version:
+        //   let delta_e = nabla_h * (self.time_step * EPS0_INV);
+        let sigma = self.mat.conductivity();
+        let e_prev = self.in_e;
+        let delta_e = (nabla_h - e_prev.elem_mul(sigma)).elem_div(
+            sigma*self.time_step + Vec3::uniform(twice_eps0)
+        )*(R::two()*self.time_step);
+        // println!("spirv-step_e delta_e: {:?}", delta_e);
+        e_prev + delta_e + self.stim_e
+    }
+}
+
+
+pub struct StepHContext<'a, R, M> {
+    pub inv_feature_size: R,
+    pub time_step: R,
+    pub stim_h: Vec3<R>,
+    pub mat: &'a M,
+    /// Input field sampled near this location
+    pub in_e: VolumeSamplePos<R>,
+    pub in_h: Vec3<R>,
+    pub in_m: Vec3<R>,
+}
+
+impl<'a, R: Real, M: Material<R>> StepHContext<'a, R, M> {
+    pub fn step_flat_view<RM, RF, WF>(
+        meta: SimMeta<R>,
+        mat: &RM,
+        stim_h: &RF,
+        e: &RF,
+        h: &mut WF,
+        m: &mut WF,
+        idx: Vec3u,
+    )
+    where
+        RM: Index<usize, Output=M> + ?Sized,
+        RF: Index<usize, Output=Vec3<R>> + ?Sized,
+        WF: Index<usize, Output=Vec3<R>> + IndexMut<usize> + ?Sized,
+    {
+        let dim = meta.dim;
+        let stim_h_matrix = DimensionedSlice::new(dim, stim_h);
+        let mat_matrix = DimensionedSlice::new(dim, mat);
+        let e_matrix = DimensionedSlice::new(dim, e);
+        let mut h_matrix = DimensionedSlice::new(dim, h);
+        let mut m_matrix = DimensionedSlice::new(dim, m);
+
+        let stim_h = stim_h_matrix[idx];
+        let mat = &mat_matrix[idx];
+        let in_e = VolumeSamplePos::from_indexable(&e_matrix, dim, idx);
+        let in_h = h_matrix[idx];
+        let in_m = m_matrix[idx];
+
+        let update_state = StepHContext {
+            inv_feature_size: meta.inv_feature_size,
+            time_step: meta.time_step,
+            stim_h,
+            mat,
+            in_e,
+            in_h,
+            in_m,
+        };
+        let (new_h, new_m) = update_state.step_h();
+        h_matrix[idx] = new_h;
+        m_matrix[idx] = new_m;
+    }
+
+    pub fn step_h(self) -> (Vec3<R>, Vec3<R>) {
+        let mu0 = R::mu0();
+        let mu0_inv = R::mu0_inv();
+        let deltas = self.in_e.delta_e();
+        // println!("spirv-step_h delta_e_struct: {:?}", deltas);
+        // \nabla x E
+        let nabla_e = deltas.nabla() * self.inv_feature_size;
+        // println!("spirv-step_h nabla_e: {:?}", nabla_e);
+        let delta_b = nabla_e * (-self.time_step);
+
+        // Relation between these is: B = mu0*(H + M)
+        let old_h = self.in_h;
+        let old_m = self.in_m;
+        let old_b = (old_h + old_m) * mu0;
+
+        let new_b = old_b + delta_b + self.stim_h * mu0;
+        let mat = self.mat;
+        let new_m = mat.move_b_vec(old_m, new_b);
+        let new_h = new_b * mu0_inv - new_m;
+        // println!("spirv-step_h delta_h: {:?}", delta_h);
+        (new_h, new_m)
+    }
+}

crates/cross/src/vec.rs

@@ -124,13 +124,17 @@ impl<R: Real> Vec2<R> {
         self.mag_sq().sqrt()
     }
 
-    pub fn with_mag(&self, new_mag: R) -> Self {
+    pub fn with_mag(&self, new_mag: R) -> Option<Self> {
         if new_mag.is_zero() {
             // avoid div-by-zero if self.mag() == 0 and new_mag == 0
-            Vec2::new(R::zero(), R::zero())
+            Some(Self::zero())
         } else {
-            let scale = new_mag / self.mag();
-            *self * scale
+            let old_mag = self.mag();
+            if old_mag.is_zero() {
+                None
+            } else {
+                Some(*self * (new_mag / old_mag))
+            }
         }
     }
 
@@ -329,13 +333,19 @@ impl<R: Real> Vec3<R> {
         Self::new(self.x().exp(), self.y().exp(), self.z().exp())
     }
 
-    pub fn with_mag(&self, new_mag: R) -> Self {
+    /// the only condition upon which this returns `None` is if the current magnitude is zero
+    /// and the new magnitude and NON-zero.
+    pub fn with_mag(&self, new_mag: R) -> Option<Self> {
         if new_mag.is_zero() {
             // avoid div-by-zero if self.mag() == 0 and new_mag == 0
-            Self::zero()
+            Some(Self::zero())
         } else {
-            let scale = new_mag / self.mag();
-            *self * scale
+            let old_mag = self.mag();
+            if old_mag.is_zero() {
+                None
+            } else {
+                Some(*self * (new_mag / old_mag))
+            }
         }
     }
 
@@ -344,7 +354,7 @@ impl<R: Real> Vec3<R> {
         if *self == Self::zero() {
             *self
         } else {
-            self.with_mag(R::one())
+            self.with_mag(R::one()).unwrap()
         }
     }
     pub fn round(&self) -> Self {

crates/cross/src/vecu.rs

@@ -51,6 +51,10 @@ impl Vec3u {
     pub fn product_sum(&self) -> u32 {
         self.x * self.y * self.z
     }
+
+    pub fn product_sum_usize(&self) -> usize {
+        self.x as usize * self.y as usize * self.z as usize
+    }
 }
 
 impl From<(u32, u32, u32)> for Vec3u {
@@ -59,6 +63,12 @@ impl From<(u32, u32, u32)> for Vec3u {
     }
 }
 
+impl Into<(u32, u32, u32)> for Vec3u {
+    fn into(self) -> (u32, u32, u32) {
+        (self.x, self.y, self.z)
+    }
+}
+
 impl<R: Real> From<Vec3<R>> for Vec3u {
     fn from(v: Vec3<R>) -> Self {
         Self::new(v.x().to_f64() as _, v.y().to_f64() as _, v.z().to_f64() as _)

crates/post/src/bin/csv.rs

@@ -1,4 +1,5 @@
-// use crate::{Loader, LoaderCache};
+//! extracts Measurements from rendered .bc files and dumps them into a CSV
+
 use coremem_post::{Loader, LoaderCache};
 use std::path::PathBuf;
 use structopt::StructOpt;
@@ -15,17 +16,13 @@ fn main() {
     
     let mut frame = cache.load_first();
     for meas in frame.measurements() {
-        for key in meas.key_value(&**frame).keys() {
-            print!("\"{}\",", key);
-        }
+        print!("\"{}\",", meas.name());
     }
     println!("");
 
     loop {
         for meas in frame.measurements() {
-            for value in meas.key_value(&**frame).values() {
-                print!("\"{}\",", value);
-            }
+            print!("\"{}\",", meas.machine_readable().replace(",", "\\,"));
         }
         println!("");
 

crates/post/src/bin/decimate.rs

@@ -1,4 +1,5 @@
-// use crate::Loader;
+//! "decimates" the binary rendered form of a simulation
+//! by deleting all but every N rendered files
 use coremem_post::Loader;
 use std::fs;
 use std::path::PathBuf;

crates/post/src/bin/viewer.rs

@@ -1,8 +1,15 @@
+//! interactive CLI viewer for .bc files.
+//! navigate through the simulation over time and space,
+//! and toggle views to see more detail over material, electric, or magnetic changes.
+
 use coremem_post::{Loader, Viewer};
 use std::path::PathBuf;
+use std::io::Write as _;
 use std::time::Duration;
 use structopt::StructOpt;
+use crossterm::{cursor, QueueableCommand as _};
 use crossterm::event::{Event, KeyCode};
+use crossterm::style::{style, PrintStyledContent};
 
 #[derive(Debug, StructOpt)]
 struct Opt {
@@ -33,6 +40,12 @@ fn event_loop(mut viewer: Viewer) {
             }
         }
         viewer.navigate(time_steps, z_steps);
+
+        let mut stdout = std::io::stdout();
+        stdout.queue(cursor::MoveToColumn(30)).unwrap();
+        stdout.queue(PrintStyledContent(style("wasd=appearance; arrows,PgUp/PgDown=navigate; q=quit"))).unwrap();
+        stdout.flush().unwrap();
+
         let _ = crossterm::event::poll(Duration::from_millis(33)).unwrap();
     }
 }

crates/post/src/lib.rs

@@ -1,15 +1,14 @@
 //! Post-processing tools
-use coremem::meas::AbstractMeasurement;
+use coremem::meas;
 use coremem::render::{ColorTermRenderer, Renderer as _, RenderConfig, SerializedFrame};
-use coremem::sim::{SimState, StaticSim};
+use coremem::sim::{AbstractSim, GenericSim};
 
 use itertools::Itertools as _;
 use lru::LruCache;
 use rayon::{ThreadPool, ThreadPoolBuilder};
 use std::collections::HashSet;
 use std::fs::{DirEntry, File, read_dir};
-use std::io::{BufReader, Seek as _, SeekFrom};
-use std::ops::Deref;
+use std::io::BufReader;
 use std::path::{Path, PathBuf};
 use std::sync::Arc;
 use std::sync::mpsc::{self, Receiver, Sender};
@@ -33,25 +32,21 @@ pub type Result<T> = std::result::Result<T, Error>;
 
 pub struct Frame {
     path: PathBuf,
-    data: SerializedFrame<StaticSim>,
+    data: SerializedFrame<GenericSim<f32>>,
 }
 
 impl Frame {
-    pub fn measurements(&self) -> &[Box<dyn AbstractMeasurement>] {
+    pub fn measurements(&self) -> &[meas::Measurement] {
         &*self.data.measurements
     }
+    pub fn sim(&self) -> &GenericSim<f32> {
+        &self.data.state
+    }
     pub fn path(&self) -> &Path {
         &*self.path
     }
 }
 
-impl Deref for Frame {
-    type Target = StaticSim;
-    fn deref(&self) -> &Self::Target {
-        &self.data.state
-    }
-}
-
 #[derive(Default)]
 pub struct Loader {
     dir: PathBuf,
@@ -105,20 +100,19 @@ impl Loader {
 
     fn load(&self, path: &Path) -> Result<Frame> {
         let mut reader = BufReader::new(File::open(path).unwrap());
-        // Try to deserialize a couple different types of likely sims.
-        // TODO: would be good to drop a marker in the file to make sure we don't
-        // 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<f32>> =
-                bincode::deserialize_from(&mut reader)?;
-            Ok(SerializedFrame::to_static(data))
-        }).or_else(|_| -> Result<_> {
-            reader.seek(SeekFrom::Start(0)).unwrap();
-            let data: SerializedFrame<SimState<f64>> =
-                bincode::deserialize_from(reader)?;
-            Ok(SerializedFrame::to_static(data))
-        })?;
+        // let data = bincode::deserialize_from(&mut reader).or_else(|_| -> Result<_> {
+        //     reader.seek(SeekFrom::Start(0)).unwrap();
+        //     let data: SerializedFrame<SimState<f32>> =
+        //         bincode::deserialize_from(&mut reader)?;
+        //     Ok(SerializedFrame::to_static(data))
+        // }).or_else(|_| -> Result<_> {
+        //     reader.seek(SeekFrom::Start(0)).unwrap();
+        //     let data: SerializedFrame<SimState<f64>> =
+        //         bincode::deserialize_from(reader)?;
+        //     Ok(SerializedFrame::to_static(data))
+        // })?;
+        // TODO: try to decode a few common sim types (as above) if this fails?
+        let data = bincode::deserialize_from(&mut reader)?;
         Ok(Frame {
             path: path.into(),
             data
@@ -253,7 +247,7 @@ impl Viewer {
         let mut cache = LoaderCache::new(loader, 6, 6);
         let viewing = cache.load_first();
         Self {
-            z: viewing.depth() / 2,
+            z: viewing.sim().depth() / 2,
             viewing,
             cache,
             renderer: Default::default(),
@@ -263,7 +257,7 @@ impl Viewer {
     }
     pub fn navigate(&mut self, time_steps: isize, z_steps: i32) {
         let new_z = (self.z as i32).saturating_add(z_steps);
-        let new_z = new_z.max(0).min(self.viewing.depth() as i32 - 1) as u32;
+        let new_z = new_z.max(0).min(self.viewing.sim().depth() as i32 - 1) as u32;
         if time_steps == 0 && new_z == self.z && self.render_config == self.last_config {
             return;
         }
@@ -275,7 +269,10 @@ impl Viewer {
     }
     pub fn render(&self) {
         self.renderer.render_z_slice(
-            &**self.viewing, self.z, &self.viewing.data.measurements, self.render_config
+            self.viewing.sim(),
+            self.z,
+            &*meas::as_dyn_measurements(self.viewing.measurements()),
+            self.render_config,
         );
     }
     pub fn render_config(&mut self) -> &mut RenderConfig {

crates/spirv_backend/Cargo.toml

@@ -8,4 +8,4 @@ crate-type = ["dylib", "lib"]
 
 [dependencies]
 spirv-std = { git = "https://github.com/EmbarkStudios/rust-gpu", features = ["glam"] }  # MIT or Apache 2.0
-coremem_types = { path = "../types" }
+coremem_cross = { path = "../cross" }

crates/spirv_backend/src/adapt.rs

@@ -0,0 +1,53 @@
+use spirv_std::RuntimeArray;
+
+use coremem_cross::mat::Material;
+use coremem_cross::real::Real;
+use coremem_cross::step::{SimMeta, StepEContext, StepHContext};
+use coremem_cross::vec::{Vec3, Vec3u};
+
+use crate::support::SizedArray;
+
+pub(crate) fn step_h<R: Real, M: Material<R>>(
+    idx: Vec3u,
+    meta: &SimMeta<R>,
+    stimulus_h: &RuntimeArray<Vec3<R>>,
+    material: &RuntimeArray<M>,
+    e: &RuntimeArray<Vec3<R>>,
+    h: &mut RuntimeArray<Vec3<R>>,
+    m: &mut RuntimeArray<Vec3<R>>,
+) {
+    if idx.x() < meta.dim.x() && idx.y() < meta.dim.y() && idx.z() < meta.dim.z() {
+        let dim = meta.dim;
+        let len = dim.product_sum_usize();
+
+        let stim_h_array = unsafe { SizedArray::new(stimulus_h, len) };
+        let mat_array = unsafe { SizedArray::new(material, len) };
+        let e_array = unsafe { SizedArray::new(e, len) };
+        let mut h_array = unsafe { SizedArray::new(h, len) };
+        let mut m_array = unsafe { SizedArray::new(m, len) };
+
+        StepHContext::step_flat_view(*meta, &mat_array, &stim_h_array, &e_array, &mut h_array, &mut m_array, idx);
+    }
+}
+
+pub(crate) fn step_e<R: Real, M: Material<R>>(
+    idx: Vec3u,
+    meta: &SimMeta<R>,
+    stimulus_e: &RuntimeArray<Vec3<R>>,
+    material: &RuntimeArray<M>,
+    e: &mut RuntimeArray<Vec3<R>>,
+    h: &RuntimeArray<Vec3<R>>,
+) {
+    if idx.x() < meta.dim.x() && idx.y() < meta.dim.y() && idx.z() < meta.dim.z() {
+        let dim = meta.dim;
+        let len = dim.product_sum_usize();
+
+        let stim_e_array = unsafe { SizedArray::new(stimulus_e, len) };
+        let mat_array = unsafe { SizedArray::new(material, len) };
+        let mut e_array = unsafe { SizedArray::new(e, len) };
+        let h_array = unsafe { SizedArray::new(h, len) };
+
+        StepEContext::step_flat_view(*meta, &mat_array, &stim_e_array, &mut e_array, &h_array, idx);
+    }
+}
+

crates/spirv_backend/src/lib.rs

@@ -8,115 +8,79 @@
 
 extern crate spirv_std;
 
-pub use spirv_std::glam;
+use spirv_std::{glam, RuntimeArray};
 #[cfg(not(target_arch = "spirv"))]
 use spirv_std::macros::spirv;
 
-pub mod mat;
-pub mod sim;
-pub mod support;
+mod adapt;
+mod support;
 
-pub use sim::{SerializedSimMeta, SerializedStepE, SerializedStepH};
-pub use support::{Optional, UnsizedArray};
+use coremem_cross::mat::{Ferroxcube3R1MH, FullyGenericMaterial, IsoConductorOr};
+use coremem_cross::real::R32;
+use coremem_cross::step::SimMeta;
+use coremem_cross::vec::{Vec3, Vec3u};
 
-use mat::{IsoConductorOr, FullyGenericMaterial};
-
-use coremem_types::mat::{Ferroxcube3R1MH, Material};
-use coremem_types::vec::{Vec3, Vec3u};
-
-type Iso3R1 = IsoConductorOr<Ferroxcube3R1MH>;
+type Iso3R1<R> = IsoConductorOr<R, Ferroxcube3R1MH>;
 
 fn glam_vec_to_internal(v: glam::UVec3) -> Vec3u {
     Vec3u::new(v.x, v.y, v.z)
 }
 
-fn step_h<M: Material<f32>>(
-    id: Vec3u,
-    meta: &SerializedSimMeta,
-    stimulus_h: &UnsizedArray<Vec3<f32>>,
-    material: &UnsizedArray<M>,
-    e: &UnsizedArray<Vec3<f32>>,
-    h: &mut UnsizedArray<Vec3<f32>>,
-    m: &mut UnsizedArray<Vec3<f32>>,
-) {
-    if id.x() < meta.dim.x() && id.y() < meta.dim.y() && id.z() < meta.dim.z() {
-        let sim_state = SerializedStepH::new(meta, stimulus_h, material, e, h, m);
-        let update_state = sim_state.index(id);
-        update_state.step_h();
-    }
+mod private {
+    pub trait Sealed {}
 }
 
-fn step_e<M: Material<f32>>(
-    id: Vec3u,
-    meta: &SerializedSimMeta,
-    stimulus_e: &UnsizedArray<Vec3<f32>>,
-    material: &UnsizedArray<M>,
-    e: &mut UnsizedArray<Vec3<f32>>,
-    h: &UnsizedArray<Vec3<f32>>,
-) {
-    if id.x() < meta.dim.x() && id.y() < meta.dim.y() && id.z() < meta.dim.z() {
-        let sim_state = SerializedStepE::new(meta, stimulus_e, material, e, h);
-
-        let update_state = sim_state.index(id);
-        update_state.step_e();
-    }
-}
-
-/// Return the step_h/step_e entry point names for the provided material
-pub fn entry_points<M: 'static>() -> Optional<(&'static str, &'static str)> {
-    use core::any::TypeId;
-    let mappings = [
-        (TypeId::of::<FullyGenericMaterial>(),
-            ("step_h_generic_material", "step_e_generic_material")
-        ),
-        (TypeId::of::<Iso3R1>(),
-            ("step_h_iso_3r1", "step_e_iso_3r1")
-        ),
-    ];
-
-    for (id, names) in mappings {
-        if id == TypeId::of::<M>() {
-            return Optional::some(names);
-        }
-    }
-    Optional::none()
+pub trait HasEntryPoints<R>: private::Sealed {
+    fn step_h() -> &'static str;
+    fn step_e() -> &'static str;
 }
 
 macro_rules! steps {
-    ($mat:ty, $step_h:ident, $step_e:ident) => {
+    ($flt:ty, $mat:ty, $step_h:ident, $step_e:ident) => {
+        impl private::Sealed for $mat { }
+        impl HasEntryPoints<$flt> for $mat {
+            fn step_h() -> &'static str {
+                stringify!($step_h)
+            }
+            fn step_e() -> &'static str {
+                stringify!($step_e)
+            }
+        }
         // LocalSize/numthreads
         #[spirv(compute(threads(4, 4, 4)))]
         pub fn $step_h(
             #[spirv(global_invocation_id)] id: glam::UVec3,
-            #[spirv(storage_buffer, descriptor_set = 0, binding = 0)] meta: &SerializedSimMeta,
+            #[spirv(storage_buffer, descriptor_set = 0, binding = 0)] meta: &SimMeta<$flt>,
             // XXX: delete this input?
-            #[spirv(storage_buffer, descriptor_set = 0, binding = 1)] _unused_stimulus_e: &UnsizedArray<Vec3<f32>>,
-            #[spirv(storage_buffer, descriptor_set = 0, binding = 2)] stimulus_h: &UnsizedArray<Vec3<f32>>,
-            #[spirv(storage_buffer, descriptor_set = 0, binding = 3)] material: &UnsizedArray<$mat>,
-            #[spirv(storage_buffer, descriptor_set = 0, binding = 4)] e: &UnsizedArray<Vec3<f32>>,
-            #[spirv(storage_buffer, descriptor_set = 0, binding = 5)] h: &mut UnsizedArray<Vec3<f32>>,
-            #[spirv(storage_buffer, descriptor_set = 0, binding = 6)] m: &mut UnsizedArray<Vec3<f32>>,
+            #[spirv(storage_buffer, descriptor_set = 0, binding = 1)] _unused_stimulus_e: &RuntimeArray<Vec3<$flt>>,
+            #[spirv(storage_buffer, descriptor_set = 0, binding = 2)] stimulus_h: &RuntimeArray<Vec3<$flt>>,
+            #[spirv(storage_buffer, descriptor_set = 0, binding = 3)] material: &RuntimeArray<$mat>,
+            #[spirv(storage_buffer, descriptor_set = 0, binding = 4)] e: &RuntimeArray<Vec3<$flt>>,
+            #[spirv(storage_buffer, descriptor_set = 0, binding = 5)] h: &mut RuntimeArray<Vec3<$flt>>,
+            #[spirv(storage_buffer, descriptor_set = 0, binding = 6)] m: &mut RuntimeArray<Vec3<$flt>>,
         ) {
-            step_h(glam_vec_to_internal(id), meta, stimulus_h, material, e, h, m)
+            adapt::step_h(glam_vec_to_internal(id), meta, stimulus_h, material, e, h, m)
         }
 
         #[spirv(compute(threads(4, 4, 4)))]
         pub fn $step_e(
             #[spirv(global_invocation_id)] id: glam::UVec3,
-            #[spirv(storage_buffer, descriptor_set = 0, binding = 0)] meta: &SerializedSimMeta,
-            #[spirv(storage_buffer, descriptor_set = 0, binding = 1)] stimulus_e: &UnsizedArray<Vec3<f32>>,
+            #[spirv(storage_buffer, descriptor_set = 0, binding = 0)] meta: &SimMeta<$flt>,
+            #[spirv(storage_buffer, descriptor_set = 0, binding = 1)] stimulus_e: &RuntimeArray<Vec3<$flt>>,
             // XXX: delete this input?
-            #[spirv(storage_buffer, descriptor_set = 0, binding = 2)] _unused_stimulus_h: &UnsizedArray<Vec3<f32>>,
-            #[spirv(storage_buffer, descriptor_set = 0, binding = 3)] material: &UnsizedArray<$mat>,
-            #[spirv(storage_buffer, descriptor_set = 0, binding = 4)] e: &mut UnsizedArray<Vec3<f32>>,
-            #[spirv(storage_buffer, descriptor_set = 0, binding = 5)] h: &UnsizedArray<Vec3<f32>>,
+            #[spirv(storage_buffer, descriptor_set = 0, binding = 2)] _unused_stimulus_h: &RuntimeArray<Vec3<$flt>>,
+            #[spirv(storage_buffer, descriptor_set = 0, binding = 3)] material: &RuntimeArray<$mat>,
+            #[spirv(storage_buffer, descriptor_set = 0, binding = 4)] e: &mut RuntimeArray<Vec3<$flt>>,
+            #[spirv(storage_buffer, descriptor_set = 0, binding = 5)] h: &RuntimeArray<Vec3<$flt>>,
             // XXX: can/should this m input be deleted?
-            #[spirv(storage_buffer, descriptor_set = 0, binding = 6)] _unused_m: &UnsizedArray<Vec3<f32>>,
+            #[spirv(storage_buffer, descriptor_set = 0, binding = 6)] _unused_m: &RuntimeArray<Vec3<$flt>>,
         ) {
-            step_e(glam_vec_to_internal(id), meta, stimulus_e, material, e, h)
+            adapt::step_e(glam_vec_to_internal(id), meta, stimulus_e, material, e, h)
         }
     };
 }
 
-steps!(FullyGenericMaterial, step_h_generic_material, step_e_generic_material);
-steps!(Iso3R1, step_h_iso_3r1, step_e_iso_3r1);
+steps!(f32, FullyGenericMaterial<f32>, step_h_generic_material_f32, step_e_generic_material_f32);
+steps!(f32, Iso3R1<f32>, step_h_iso_3r1_f32, step_e_iso_3r1_f32);
+steps!(R32, FullyGenericMaterial<R32>, step_h_generic_material_r32, step_e_generic_material_r32);
+steps!(R32, Iso3R1<R32>, step_h_iso_3r1_r32, step_e_iso_3r1_r32);

crates/spirv_backend/src/mat.rs

@@ -1,4 +0,0 @@
-// TODO: plumb the R parameter through and remove this.
-pub type FullyGenericMaterial = coremem_types::mat::FullyGenericMaterial<f32>;
-pub type IsoConductorOr<M> = coremem_types::mat::IsoConductorOr<f32, M>;
-

crates/spirv_backend/src/sim.rs

@@ -1,338 +0,0 @@
-// use spirv_std::RuntimeArray;
-use crate::support::{
-    Array3, Array3Mut, ArrayHandle, ArrayHandleMut, Optional, UnsizedArray
-};
-use coremem_types::mat::Material;
-use coremem_types::real::{Real as _};
-use coremem_types::vec::{Vec3, Vec3u};
-
-#[derive(Copy, Clone)]
-pub struct SerializedSimMeta {
-    pub dim: Vec3u,
-    pub inv_feature_size: f32,
-    pub time_step: f32,
-    pub feature_size: f32,
-}
-
-
-/// Whatever data we received from the host in their call to step_h
-pub struct SerializedStepH<'a, M> {
-    meta: &'a SerializedSimMeta,
-    stimulus_h: &'a UnsizedArray<Vec3<f32>>,
-    material: &'a UnsizedArray<M>,
-    e: &'a UnsizedArray<Vec3<f32>>,
-    h: &'a mut UnsizedArray<Vec3<f32>>,
-    m: &'a mut UnsizedArray<Vec3<f32>>,
-}
-
-impl<'a, M> SerializedStepH<'a, M> {
-    pub fn new(
-        meta: &'a SerializedSimMeta,
-        stimulus_h: &'a UnsizedArray<Vec3<f32>>,
-        material: &'a UnsizedArray<M>,
-        e: &'a UnsizedArray<Vec3<f32>>,
-        h: &'a mut UnsizedArray<Vec3<f32>>,
-        m: &'a mut UnsizedArray<Vec3<f32>>,
-    ) -> Self {
-        Self {
-            meta, stimulus_h, material, e, h, m
-        }
-    }
-
-    pub fn index(self, idx: Vec3u) -> StepHContext<'a, M> {
-        let dim = self.meta.dim;
-        let stim_h_matrix = Array3::new(self.stimulus_h, dim);
-        let mat_matrix = Array3::new(self.material, dim);
-        let e = Array3::new(self.e, dim);
-        let h = Array3Mut::new(self.h, dim);
-        let m = Array3Mut::new(self.m, dim);
-
-        let in_e = VolumeSamplePos {
-            mid: e.get(idx).unwrap(),
-            xp1: e.get(idx + Vec3u::unit_x()),
-            yp1: e.get(idx + Vec3u::unit_y()),
-            zp1: e.get(idx + Vec3u::unit_z()),
-        };
-        let out_h = h.into_mut_handle(idx);
-        let out_m = m.into_mut_handle(idx);
-
-        let mat = mat_matrix.into_handle(idx);
-
-        StepHContext {
-            inv_feature_size: self.meta.inv_feature_size,
-            time_step: self.meta.time_step,
-            stim_h: stim_h_matrix.get(idx).unwrap(),
-            mat,
-            in_e,
-            out_h,
-            out_m,
-        }
-    }
-}
-
-/// Whatever data we received from the host in their call to step_e
-pub struct SerializedStepE<'a, M> {
-    meta: &'a SerializedSimMeta,
-    stimulus_e: &'a UnsizedArray<Vec3<f32>>,
-    material: &'a UnsizedArray<M>,
-    e: &'a mut UnsizedArray<Vec3<f32>>,
-    h: &'a UnsizedArray<Vec3<f32>>,
-}
-
-impl<'a, M> SerializedStepE<'a, M> {
-    pub fn new(
-        meta: &'a SerializedSimMeta,
-        stimulus_e: &'a UnsizedArray<Vec3<f32>>,
-        material: &'a UnsizedArray<M>,
-        e: &'a mut UnsizedArray<Vec3<f32>>,
-        h: &'a UnsizedArray<Vec3<f32>>,
-    ) -> Self {
-        Self {
-            meta, stimulus_e, material, e, h
-        }
-    }
-
-    pub fn index(self, idx: Vec3u) -> StepEContext<'a, M> {
-        let dim = self.meta.dim;
-        let stim_e_matrix = Array3::new(self.stimulus_e, dim);
-        let mat_matrix = Array3::new(self.material, dim);
-        let e = Array3Mut::new(self.e, dim);
-        let h = Array3::new(self.h, dim);
-
-        let xm1 = if idx.x() == 0 {
-            Optional::none()
-        } else {
-            h.get(idx - Vec3u::unit_x())
-        };
-        let ym1 = if idx.y() == 0 {
-            Optional::none()
-        } else {
-            h.get(idx - Vec3u::unit_y())
-        };
-        let zm1 = if idx.z() == 0 {
-            Optional::none()
-        } else {
-            h.get(idx - Vec3u::unit_z())
-        };
-
-        let in_h = VolumeSampleNeg {
-            mid: h.get(idx).unwrap(),
-            xm1,
-            ym1,
-            zm1,
-        };
-        let out_e = e.into_mut_handle(idx);
-
-        let mat = mat_matrix.into_handle(idx);
-
-        StepEContext {
-            inv_feature_size: self.meta.inv_feature_size,
-            time_step: self.meta.time_step,
-            stim_e: stim_e_matrix.get(idx).unwrap(),
-            mat,
-            in_h,
-            out_e,
-        }
-    }
-}
-
-/// Package the field vectors adjacent to some particular location.
-/// Particular those at negative offsets from the midpoint.
-/// This is used in step_e when looking at the H field deltas.
-#[derive(Copy, Clone)]
-struct VolumeSampleNeg {
-    mid: Vec3<f32>,
-    xm1: Optional<Vec3<f32>>,
-    ym1: Optional<Vec3<f32>>,
-    zm1: Optional<Vec3<f32>>,
-}
-
-impl VolumeSampleNeg {
-    /// Calculate the delta in H values amongst this cell and its neighbors (left/up/out)
-    fn delta_h(self) -> FieldDeltas {
-        let mid = self.mid;
-        // let (dfy_dx, dfz_dx) = self.xm1.map(|xm1| {
-        //     (mid.y() - xm1.y(), mid.z() - xm1.z())
-        // }).unwrap_or_default();
-
-        // let (dfx_dy, dfz_dy) = self.ym1.map(|ym1| {
-        //     (mid.x() - ym1.x(), mid.z() - ym1.z())
-        // }).unwrap_or_default();
-
-        // let (dfx_dz, dfy_dz) = self.zm1.map(|zm1| {
-        //     (mid.x() - zm1.x(), mid.y() - zm1.y())
-        // }).unwrap_or_default();
-
-        let (dfy_dx, dfz_dx) = if self.xm1.is_some() {
-            (mid.y() - self.xm1.unwrap().y(), mid.z() - self.xm1.unwrap().z())
-        } else {
-            (0.0, 0.0)
-        };
-
-        let (dfx_dy, dfz_dy) = if self.ym1.is_some() {
-            (mid.x() - self.ym1.unwrap().x(), mid.z() - self.ym1.unwrap().z())
-        } else {
-            (0.0, 0.0)
-        };
-
-        let (dfx_dz, dfy_dz) = if self.zm1.is_some() {
-            (mid.x() - self.zm1.unwrap().x(), mid.y() - self.zm1.unwrap().y())
-        } else {
-            (0.0, 0.0)
-        };
-        
-        FieldDeltas {
-            dfy_dx,
-            dfz_dx,
-            dfx_dy,
-            dfz_dy,
-            dfx_dz,
-            dfy_dz,
-        }
-    }
-}
-
-
-/// Package the field vectors adjacent to some particular location.
-/// Particular those at positive offsets from the midpoint.
-/// This is used in step_h when looking at the E field deltas.
-#[derive(Copy, Clone)]
-struct VolumeSamplePos {
-    mid: Vec3<f32>,
-    xp1: Optional<Vec3<f32>>,
-    yp1: Optional<Vec3<f32>>,
-    zp1: Optional<Vec3<f32>>
-}
-
-impl VolumeSamplePos {
-    /// Calculate the delta in E values amongst this cell and its neighbors (right/down/in)
-    fn delta_e(self) -> FieldDeltas {
-        let mid = self.mid;
-        // let (dfy_dx, dfz_dx) = self.xp1.map(|xp1| {
-        //     (xp1.y() - mid.y(), xp1.z() - mid.z())
-        // }).unwrap_or_default();
-
-        // let (dfx_dy, dfz_dy) = self.yp1.map(|yp1| {
-        //     (yp1.x() - mid.x(), yp1.z() - mid.z())
-        // }).unwrap_or_default();
-
-        // let (dfx_dz, dfy_dz) = self.zp1.map(|zp1| {
-        //     (zp1.x() - mid.x(), zp1.y() - mid.y())
-        // }).unwrap_or_default();
-
-        let (dfy_dx, dfz_dx) = if self.xp1.is_some() {
-            (self.xp1.unwrap().y() - mid.y(), self.xp1.unwrap().z() - mid.z())
-        } else {
-            (0.0, 0.0)
-        };
-
-        let (dfx_dy, dfz_dy) = if self.yp1.is_some() {
-            (self.yp1.unwrap().x() - mid.x(), self.yp1.unwrap().z() - mid.z())
-        } else {
-            (0.0, 0.0)
-        };
-
-        let (dfx_dz, dfy_dz) = if self.zp1.is_some() {
-            (self.zp1.unwrap().x() - mid.x(), self.zp1.unwrap().y() - mid.y())
-        } else {
-            (0.0, 0.0)
-        };
-        
-        FieldDeltas {
-            dfy_dx,
-            dfz_dx,
-            dfx_dy,
-            dfz_dy,
-            dfx_dz,
-            dfy_dz,
-        }
-    }
-}
-
-struct FieldDeltas {
-    dfy_dx: f32,
-    dfz_dx: f32,
-    dfx_dy: f32,
-    dfz_dy: f32,
-    dfx_dz: f32,
-    dfy_dz: f32,
-}
-
-impl FieldDeltas {
-    fn nabla(self) -> Vec3<f32> {
-        Vec3::new(
-            self.dfz_dy - self.dfy_dz,
-            self.dfx_dz - self.dfz_dx,
-            self.dfy_dx - self.dfx_dy,
-        )
-    }
-}
-
-pub struct StepEContext<'a, M> {
-    inv_feature_size: f32,
-    time_step: f32,
-    stim_e: Vec3<f32>,
-    mat: ArrayHandle<'a, M>,
-    /// Input field sampled near this location
-    in_h: VolumeSampleNeg,
-    /// Handle to the output field at one specific index.
-    out_e: ArrayHandleMut<'a, Vec3<f32>>,
-}
-
-impl<'a, M: Material<f32>> StepEContext<'a, M> {
-    pub fn step_e(mut self) {
-        let twice_eps0 = f32::twice_eps0();
-        let deltas = self.in_h.delta_h();
-        // \nabla x H
-        let nabla_h = deltas.nabla() * self.inv_feature_size;
-        // $\nabla x H = \epsilon_0 dE/dt + \sigma E$
-        // no-conductivity version:
-        //   let delta_e = nabla_h * (self.time_step * EPS0_INV);
-        let sigma = self.mat.get_ref().conductivity();
-        let e_prev = self.out_e.get();
-        let delta_e = (nabla_h - e_prev.elem_mul(sigma)).elem_div(
-            sigma*self.time_step + Vec3::uniform(twice_eps0)
-        )*(2.0*self.time_step);
-        // println!("spirv-step_e delta_e: {:?}", delta_e);
-        self.out_e.write(e_prev + delta_e + self.stim_e);
-    }
-}
-
-
-pub struct StepHContext<'a, M> {
-    inv_feature_size: f32,
-    time_step: f32,
-    stim_h: Vec3<f32>,
-    mat: ArrayHandle<'a, M>,
-    /// Input field sampled near this location
-    in_e: VolumeSamplePos,
-    /// Handle to the output field at one specific index.
-    out_h: ArrayHandleMut<'a, Vec3<f32>>,
-    out_m: ArrayHandleMut<'a, Vec3<f32>>,
-}
-
-impl<'a, M: Material<f32>> StepHContext<'a, M> {
-    pub fn step_h(mut self) {
-        let mu0 = f32::mu0();
-        let mu0_inv = f32::mu0_inv();
-        let deltas = self.in_e.delta_e();
-        // println!("spirv-step_h delta_e_struct: {:?}", deltas);
-        // \nabla x E
-        let nabla_e = deltas.nabla() * self.inv_feature_size;
-        // println!("spirv-step_h nabla_e: {:?}", nabla_e);
-        let delta_b = nabla_e * (-self.time_step);
-
-        // Relation between these is: B = mu0*(H + M)
-        let old_h = self.out_h.get();
-        let old_m = self.out_m.get();
-        let old_b = (old_h + old_m) * mu0;
-
-        let new_b = old_b + delta_b + self.stim_h * mu0;
-        let mat = self.mat.get_ref();
-        let new_m = mat.move_b_vec(old_m, new_b);
-        let new_h = new_b * mu0_inv - new_m;
-        // println!("spirv-step_h delta_h: {:?}", delta_h);
-        self.out_h.write(new_h);
-        self.out_m.write(new_m);
-    }
-}

crates/spirv_backend/src/support.rs

@@ -1,284 +1,51 @@
-use coremem_types::vec::Vec3u;
+use core::ops::{Index, IndexMut};
 
-/// This is a spirv-compatible option type.
-/// The native rust Option type produces invalid spirv due to its enum nature; this custom option
-/// type creates code which will actually compile.
-#[derive(Copy, Clone, PartialEq)]
-pub struct Optional<T> {
-    // XXX: not a bool, because: "entrypoint parameter cannot contain a boolean"
-    present: u8,
-    data: T,
+use spirv_std::RuntimeArray;
+
+/// this is intended to wrap an unsized array with a length and provide safe indexing operations
+/// into it. it behaves quite similar to a native rust slice, and ideally we'd just use that but
+/// spirv support for slices is poor.
+pub struct SizedArray<T> {
+    items: T,
+    len: usize,
 }
 
-impl<T> Optional<T> {
-    pub fn some(data: T) -> Self {
-        Self {
-            present: 1,
-            data,
-        }
-    }
-
-    pub fn explicit_none(data: T) -> Self {
-        Self {
-            present: 0,
-            data,
-        }
-    }
-
-    pub fn is_some(self) -> bool {
-        self.present != 0
-    }
-
-    pub fn unwrap(self) -> T {
-        assert!(self.present != 0);
-        self.data
-    }
-
-    pub fn map<U: Default, F: FnOnce(T) -> U>(self, f: F) -> Optional<U> {
-        self.and_then(|inner| Optional::some(f(inner)))
-    }
-
-    pub fn and_then<U: Default, F: FnOnce(T) -> Optional<U>>(self, f: F) -> Optional<U> {
-        if self.present != 0 {
-            f(self.unwrap())
-        } else {
-            Optional::none()
-        }
-    }
-
-    pub fn unwrap_or(self, default: T) -> T {
-        if self.present != 0 {
-            self.unwrap()
-        } else {
-            default
-        }
+impl<T> SizedArray<T> {
+    /// construct the slice from some other collection and the number of items in said collection.
+    /// safety: caller must validate that it's safe to index all `len` elements in the underlying
+    /// collection.
+    pub unsafe fn new(items: T, len: usize) -> Self {
+        Self { items, len }
     }
 }
 
-impl<T: Default> Optional<T> {
-    pub fn none() -> Self {
-        Self::explicit_none(Default::default())
-    }
+impl<'a, T> Index<usize> for SizedArray<&'a RuntimeArray<T>> {
+    type Output=T;
 
-    pub fn unwrap_or_default(self) -> T {
-        self.unwrap_or(Default::default())
-    }
-}
-
-impl<T: Default> Default for Optional<T> {
-    fn default() -> Self {
-        Self::none()
-    }
-}
-
-impl<T0: Default, T1: Default> Optional<(T0, T1)> {
-    pub fn flatten((f0, f1): (Optional<T0>, Optional<T1>)) -> Self {
-        if f0.present != 0 && f1.present != 0 {
-            Optional::some((f0.unwrap(), f1.unwrap()))
-        } else {
-            Optional::none()
-        }
-    }
-}
-
-/// This struct allows doing things like *(ptr + offset) = value.
-/// Such code wouldn't ordinarily compile with the spirv target because of
-/// unsupported pointer math. RuntimeArray exists to overcome this, however
-/// it emits invalid code for non-primitive types. Hence, this hack.
-// XXX: maximum bytes an array may occupy is 0x7fff_ffff
-// We don't know the element size, so assume it to be <= 44 bytes
-// pub const MAX_UNSIZED_ARRAY: usize = 0x2e7_ffff;
-pub const MAX_UNSIZED_ARRAY: usize = 0x1ff_ffff;
-pub struct UnsizedArray<T>(pub [T; MAX_UNSIZED_ARRAY]);
-
-impl<T: Copy> UnsizedArray<T> {
-    pub unsafe fn index(&self, index: usize) -> T {
-        // *self.0.index(index)
-        self.0[index]
-        // *self.0.get_unchecked(index)
-        // &self.0
-        //asm! {
-        //    "%result = OpAccessChain _ {arr} {index}",
-        //    "OpReturnValue %result",
-        //    "%unused = OpLabel",
-        //    arr = in(reg) self,
-        //    index = in(reg) index,
-        //}
-    }
-
-    pub unsafe fn write(&mut self, index: usize, value: T) {
-        self.0[index] = value;
-    }
-}
-
-impl<T> UnsizedArray<T> {
-    pub unsafe fn index_ref(&self, index: usize) -> &T {
-        &self.0[index]
-    }
-
-    pub unsafe fn get_handle<'a>(&'a self, index: usize) -> ArrayHandle<'a, T> {
-        ArrayHandle {
-            base: self,
-            index,
-        }
-    }
-
-    pub unsafe fn get_handle_mut<'a>(&'a mut self, index: usize) -> ArrayHandleMut<'a, T> {
-        ArrayHandleMut {
-            base: self,
-            index,
-        }
-    }
-}
-
-pub struct ArrayHandle<'a, T> {
-    base: &'a UnsizedArray<T>,
-    index: usize,
-}
-
-impl<'a, T: Copy> ArrayHandle<'a, T> {
-    pub fn get(&self) -> T {
+    fn index(&self, idx: usize) -> &Self::Output {
+        debug_assert!(idx < self.len);
         unsafe {
-            self.base.index(self.index)
+            self.items.index(idx)
         }
     }
-    pub fn get_into(&self, out: &mut T) {
-        core::clone::Clone::clone_from(out, self.get_ref());
-    }
 }
-impl<'a, T> ArrayHandle<'a, T> {
-    pub fn get_ref(&self) -> &T {
+
+impl<'a, T> Index<usize> for SizedArray<&'a mut RuntimeArray<T>> {
+    type Output=T;
+
+    fn index(&self, idx: usize) -> &Self::Output {
+        debug_assert!(idx < self.len);
         unsafe {
-            self.base.index_ref(self.index)
+            self.items.index(idx)
         }
     }
 }
-
-
-pub struct ArrayHandleMut<'a, T> {
-    base: &'a mut UnsizedArray<T>,
-    index: usize,
-}
-
-impl<'a, T: Copy> ArrayHandleMut<'a, T> {
-    pub fn write(&mut self, value: T) {
+impl<'a, T> IndexMut<usize> for SizedArray<&'a mut RuntimeArray<T>> {
+    fn index_mut(&mut self, idx: usize) -> &mut Self::Output {
+        debug_assert!(idx < self.len);
         unsafe {
-            self.base.write(self.index, value);
-        }
-    }
-    pub fn get(&self) -> T {
-        unsafe {
-            self.base.index(self.index)
+            self.items.index_mut(idx)
         }
     }
 }
 
-/// 3d dynamically-sized array backed by a borrowed buffer.
-#[derive(Clone, Copy)]
-pub struct Array3<'a, T> {
-    data: &'a UnsizedArray<T>,
-    dim: Vec3u,
-}
-
-fn index(loc: Vec3u, dim: Vec3u) -> usize {
-    ((loc.z()*dim.y() + loc.y())*dim.x() + loc.x()) as usize
-}
-
-fn checked_index(idx: Vec3u, dim: Vec3u) -> Optional<usize> {
-    if idx.x() < dim.x() && idx.y() < dim.y() && idx.z() < dim.z() {
-        let flat_idx = index(idx, dim);
-        Optional::some(flat_idx)
-    } else {
-        Optional::none()
-    }
-}
-
-impl<'a, T> Array3<'a, T> {
-    pub fn new(data: &'a UnsizedArray<T>, dim: Vec3u) -> Self {
-        Self {
-            data,
-            dim,
-        }
-    }
-
-    pub fn index(self, idx: Vec3u) -> Optional<usize> {
-        checked_index(idx, self.dim)
-    }
-
-    pub fn into_handle(self, idx: Vec3u) -> ArrayHandle<'a, T> {
-        let idx = checked_index(idx, self.dim).unwrap();
-        unsafe {
-            self.data.get_handle(idx)
-        }
-    }
-}
-
-impl<'a, T: Copy + Default> Array3<'a, T> {
-    pub fn get(&self, idx: Vec3u) -> Optional<T> {
-        let idx = self.index(idx);
-        if idx.is_some() {
-            Optional::some(unsafe {
-                self.data.index(idx.unwrap())
-            })
-        } else {
-            Optional::none()
-        }
-    }
-}
-
-/// 3d dynamically-sized array backed by a mutably borrowed buffer.
-pub struct Array3Mut<'a, T> {
-    data: &'a mut UnsizedArray<T>,
-    dim: Vec3u,
-}
-
-impl<'a, T> Array3Mut<'a, T> {
-    pub fn new(data: &'a mut UnsizedArray<T>, dim: Vec3u) -> Self {
-        Self {
-            data,
-            dim,
-        }
-    }
-
-    // fn as_array3<'b>(&'b self) -> Array3<'b, T> {
-    //     Array3::new(self.data, self.dim)
-    // }
-
-    pub fn index(&self, idx: Vec3u) -> Optional<usize> {
-        if idx.x() < self.dim.x() && idx.y() < self.dim.y() && idx.z() < self.dim.z() {
-            let flat_idx = index(idx, self.dim);
-            Optional::some(flat_idx)
-        } else {
-            Optional::none()
-        }
-    }
-}
-
-impl<'a, T: Copy> Array3Mut<'a, T> {
-    pub fn into_mut_handle(self, idx: Vec3u) -> ArrayHandleMut<'a, T> {
-        let idx = self.index(idx).unwrap();
-        unsafe {
-            self.data.get_handle_mut(idx)
-        }
-    }
-}
-
-
-#[cfg(test)]
-mod test {
-    use super::*;
-    #[test]
-    fn test_index() {
-        let dim = Vec3u::new(2, 3, 7);
-        assert_eq!(index(Vec3u::new(0, 0, 0), dim), 0);
-        assert_eq!(index(Vec3u::new(1, 0, 0), dim), 1);
-        assert_eq!(index(Vec3u::new(0, 1, 0), dim), 2);
-        assert_eq!(index(Vec3u::new(1, 1, 0), dim), 3);
-        assert_eq!(index(Vec3u::new(0, 2, 0), dim), 4);
-        assert_eq!(index(Vec3u::new(0, 0, 1), dim), 6);
-        assert_eq!(index(Vec3u::new(1, 0, 1), dim), 7);
-        assert_eq!(index(Vec3u::new(0, 1, 1), dim), 8);
-        assert_eq!(index(Vec3u::new(1, 2, 1), dim), 11);
-        assert_eq!(index(Vec3u::new(1, 2, 2), dim), 17);
-    }
-}

crates/types/src/compound/peano.rs

@@ -1,46 +0,0 @@
-#[derive(Copy, Clone, Default, PartialEq)]
-pub struct PNext<P>(P);
-#[derive(Copy, Clone, Default, PartialEq)]
-pub struct P0;
-
-pub type P1 = PNext<P0>;
-pub type P2 = PNext<P1>;
-pub type P3 = PNext<P2>;
-pub type P4 = PNext<P3>;
-pub type P5 = PNext<P4>;
-pub type P6 = PNext<P5>;
-pub type P7 = PNext<P6>;
-
-pub trait Peano: Copy + Clone + Default + PartialEq {
-    type Next: PeanoNonZero;
-    type PrevOrZero: Peano;
-    /// always set to ()
-    /// this exists to allow Peano numbers to be used as struct parameters without PhantomData
-    type Unit: Copy + Clone + Default + PartialEq;
-    const VALUE: u32;
-}
-
-pub trait PeanoNonZero: Peano {
-    type Prev: Peano;
-}
-
-impl Peano for P0 {
-    type Next = P1;
-    type PrevOrZero = P0;
-    type Unit = ();
-    const VALUE: u32 = 0;
-}
-impl<P: Peano> Peano for PNext<P> {
-    type Next = PNext<PNext<P>>;
-    type PrevOrZero = P;
-    type Unit = ();
-    const VALUE: u32 = 1 + P::VALUE;
-}
-impl<P: Peano> PeanoNonZero for PNext<P> {
-    type Prev = P;
-}
-
-// A: LessThan<B> is satisfied only if A is strictly less than B.
-pub trait LessThan<P: Peano> { }
-impl<P: Peano> LessThan<PNext<P>> for P { }
-

crates/types/src/lib.rs

@@ -1,8 +0,0 @@
-#![no_std]
-#![feature(core_intrinsics)]
-
-mod compound;
-pub mod mat;
-pub mod real;
-pub mod vec;
-mod vecu;