diff --git a/examples/minimal_torus.rs b/examples/minimal_torus.rs
new file mode 100644
index 0000000..78f082a
--- /dev/null
+++ b/examples/minimal_torus.rs
@@ -0,0 +1,99 @@
+use coremem::{Driver, Flt, mat, meas};
+use coremem::geom::{Index, Meters, Torus};
+use coremem::stim::{CurlStimulus, Sinusoid1, TimeVarying1 as _};
+
+fn main() {
+    coremem::init_logging();
+    let feat_size = 10e-6; // feature size
+    let duration = 3e-11;
+    let width = 2200e-6;
+    let depth = 1800e-6;
+    let buffer = 200e-6;
+    let ferro_major = 320e-6;
+    let ferro_minor = 60e-6;
+    let wire_minor = 40e-6;
+    let wire_major = 160e-6;
+    let peak_current = 5e11;
+    let current_duration = 1.0e-11; // half-wavelength of the sine wave
+    let current_break = 0.5e-11; // time between 'set' pulse and 'clear' pulse
+    let conductivity = 1.0e9;
+
+    let from_m = |m: Flt| (m/feat_size).round() as u32;
+    let m_to_um = |m: Flt| (m * 1e6).round() as u32;
+    let half_width = width * 0.5;
+    let half_depth = depth * 0.5;
+
+    let width_px = from_m(width);
+    let depth_px = from_m(depth);
+    let size_px = Index((width_px, width_px, depth_px).into());
+    let mut driver: Driver<mat::GenericMaterial> = Driver::new(size_px, feat_size);
+    driver.set_steps_per_frame(10);
+    driver.set_steps_per_stim(1);
+    let base = "minimal_torus-2";
+    let ferro_region = Torus::new_xy(Meters((half_width, half_width, half_depth).into()), ferro_major, ferro_minor);
+    let drive_region = Torus::new_xz(Meters((half_width - ferro_major, half_width, half_depth).into()), wire_major, wire_minor);
+    let sense_region = Torus::new_xz(Meters((half_width + ferro_major, half_width, half_depth).into()), wire_major, wire_minor);
+
+    driver.fill_region(&ferro_region, mat::db::minimal_square_ferrite().into());
+    driver.fill_region(&drive_region, mat::db::conductor(conductivity).into());
+    driver.fill_region(&sense_region, mat::db::conductor(conductivity).into());
+
+    let boundary_xy = half_width - ferro_major - ferro_minor - buffer;
+    let boundary_z = half_depth - wire_major - wire_minor - buffer;
+    println!("boundary: {}um; {}um", m_to_um(boundary_xy), m_to_um(boundary_z));
+    driver.add_upml_boundary(Meters((boundary_xy, boundary_xy, boundary_z).into()));
+
+    // J=\sigma E
+    // dJ/dt = \sigma dE/dT
+    // dE/dt = dJ/dt / \sigma
+    // dE/dt = dI/dt / (A*\sigma)
+    // if I = k*sin(w t) then dE/dt = k*w sin(w t) / (A*\sigma)
+    // i.e. dE/dt is proportional to I/(A*\sigma), multiplied by w (or, divided by wavelength)
+    let peak_stim = peak_current/current_duration / (drive_region.cross_section() * conductivity);
+    let pos_wave = Sinusoid1::from_wavelength(peak_stim, current_duration * 2.0)
+        .half_cycle();
+
+    let neg_wave = Sinusoid1::from_wavelength(-peak_stim, current_duration * 2.0)
+        .half_cycle()
+        .shifted(current_duration + current_break);
+
+    driver.add_stimulus(CurlStimulus::new(
+        drive_region.clone(),
+        pos_wave,
+        drive_region.center(),
+        drive_region.axis()
+    ));
+    driver.add_stimulus(CurlStimulus::new(
+        drive_region.clone(),
+        neg_wave,
+        drive_region.center(),
+        drive_region.axis()
+    ));
+
+    driver.add_measurement(meas::CurrentLoop::new("sense", sense_region.clone()));
+    driver.add_measurement(meas::Current::new("sense", sense_region.clone()));
+    driver.add_measurement(meas::MagneticLoop::new("mem", ferro_region.clone()));
+    driver.add_measurement(meas::Magnetization::new("mem", ferro_region.clone()));
+    driver.add_measurement(meas::MagneticFlux::new("mem", ferro_region.clone()));
+    driver.add_measurement(meas::CurrentLoop::new("drive", drive_region.clone()));
+    driver.add_measurement(meas::Current::new("drive", drive_region.clone()));
+
+    let prefix = format!("out/{}/{}-flt{}-{}-feat{}um-{}mA-{}ps--radii{}um-{}um-{}um-{}um",
+        base,
+        base,
+        std::mem::size_of::<Flt>() * 8,
+        *size_px,
+        m_to_um(feat_size),
+        (peak_current * 1e3).round() as i64,
+        (current_duration * 1e12).round() as i64,
+        m_to_um(ferro_major),
+        m_to_um(ferro_minor),
+        m_to_um(wire_major),
+        m_to_um(wire_minor),
+    );
+    let _ = std::fs::create_dir_all(&prefix);
+
+    driver.add_serializer_renderer(&*format!("{}/frame-", prefix));
+
+    driver.step_until(duration);
+}
diff --git a/src/driver.rs b/src/driver.rs
index c76e989..7867d6e 100644
--- a/src/driver.rs
+++ b/src/driver.rs
@@ -24,6 +24,7 @@ pub struct Driver<M=GenericMaterial> {
     steps_per_frame: u64,
     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>>,
@@ -42,6 +43,7 @@ impl<M: Default> Driver<M> {
             steps_per_frame: 1,
             time_spent_stepping: Default::default(),
             time_spent_on_stimuli: Default::default(),
+            time_spent_prepping_render: Default::default(),
             time_spent_blocked_on_render: Default::default(),
             time_spent_rendering: Default::default(),
             measurements: vec![
@@ -124,6 +126,7 @@ impl<M: Material + Serialize + Clone + Send + Sync + 'static> Driver<M> {
 
 impl<M: Material + Clone + Default + Send + Sync + 'static> Driver<M> {
     fn render(&mut self) {
+        let prep_start = Instant::now();
         let their_state = self.state.clone();
         let their_measurements = self.measurements.clone();
         let renderer = self.renderer.clone();
@@ -138,6 +141,7 @@ impl<M: Material + Clone + Default + Send + Sync + 'static> Driver<M> {
             *time_spent_rendering.lock().unwrap() += start_time.elapsed();
             trace!("render end");
         });
+        self.time_spent_prepping_render += prep_start.elapsed();
         let block_start = Instant::now();
         self.render_channel.1.recv().unwrap();
         self.time_spent_blocked_on_render += block_start.elapsed();
@@ -166,12 +170,13 @@ impl<M: Material + Clone + Default + Send + Sync + 'static> Driver<M> {
             let step_time = self.time_spent_stepping.as_secs_f64();
             let stim_time = self.time_spent_on_stimuli.as_secs_f64();
             let render_time = self.time_spent_rendering.lock().unwrap().as_secs_f64();
+            let render_prep_time = self.time_spent_prepping_render.as_secs_f64();
             let block_time = self.time_spent_blocked_on_render.as_secs_f64();
             let overall_time = self.start_time.elapsed().as_secs_f64();
             let fps = (self.state.step_no() as f64) / overall_time;
             let sim_time = self.state.time() as f64;
             info!(
-                "t={:.2e} frame {:06} fps: {:6.2} (sim: {:.1}s, stim: {:.1}s, render: {:.1}s, blocked: {:.1}s, other: {:.1}s)",
+                "t={:.2e} frame {:06} fps: {:6.2} (sim: {:.1}s, stim: {:.1}s, [render: {:.1}s], blocked: {:.1}s, render_prep: {:.1}s, other: {:.1}s)",
                 sim_time,
                 step,
                 fps,
@@ -179,7 +184,8 @@ impl<M: Material + Clone + Default + Send + Sync + 'static> Driver<M> {
                 stim_time,
                 render_time,
                 block_time,
-                overall_time - step_time - stim_time - render_time
+                render_prep_time,
+                overall_time - step_time - stim_time - block_time - render_prep_time,
             );
         }
     }
diff --git a/src/mat.rs b/src/mat.rs
index b00231d..25a24c4 100644
--- a/src/mat.rs
+++ b/src/mat.rs
@@ -118,6 +118,7 @@ impl MHCurve {
             geom: Polygon2d::new(full_pts)
         }
     }
+
     fn from_bh(points: &[(Flt, Flt)]) -> Self {
         let mh_points: Vec<Vec2> = points.iter().cloned().map(|(h, b)| {
             Vec2::new(h, b / consts::MU0 - h)
@@ -126,6 +127,14 @@ impl MHCurve {
         Self::new(&*mh_points)
     }
 
+    fn from_mh(points: &[(Flt, Flt)]) -> Self {
+        let mh_points: Vec<Vec2> = points.iter().cloned().map(|(h, m)| {
+            Vec2::new(h, m)
+        }).collect();
+
+        Self::new(&*mh_points)
+    }
+
     /// Return (Hmax, Mmax)
     pub fn extremes(&self) -> Vec2 {
         Vec2::new(self.geom.max_x(), self.geom.max_y())
@@ -238,11 +247,41 @@ impl PiecewiseLinearFerromagnet for Ferroxcube3R1 {
     }
 }
 
+#[derive(Default, Copy, Clone, Serialize, Deserialize)]
+pub struct MinimalSquare {
+    m: Vec3,
+}
+
+impl PiecewiseLinearFerromagnet for MinimalSquare {
+    fn curve() -> &'static MHCurve {
+        lazy_static! {
+            static ref CURVE: MHCurve = MHCurve::from_mh(&[
+                (  1.0,       0.0),
+                (  2.0, 1000000.0),
+                // Falling
+                (  0.0,  900000.0),
+            ]);
+        }
+        &*CURVE
+    }
+    fn conductivity() -> Flt {
+        1e-3
+    }
+    fn m(&self) -> Vec3 {
+        self.m
+    }
+    fn m_mut(&mut self) -> &mut Vec3 {
+        &mut self.m
+    }
+}
+
+
 #[enum_dispatch(Material)]
 #[derive(Clone, Serialize, Deserialize)]
 pub enum GenericMaterial {
     Conductor(Conductor),
     Ferroxcube3R1(Ferroxcube3R1),
+    MinimalSquare(MinimalSquare),
 }
 
 impl Default for GenericMaterial {
@@ -264,6 +303,9 @@ pub mod db {
     pub fn ferroxcube_3r1() -> Ferroxcube3R1 {
         Ferroxcube3R1::default()
     }
+    pub fn minimal_square_ferrite() -> MinimalSquare {
+        MinimalSquare::default()
+    }
 }
 
 #[cfg(test)]
diff --git a/src/render.rs b/src/render.rs
index 19cefae..d9cb081 100644
--- a/src/render.rs
+++ b/src/render.rs
@@ -102,9 +102,11 @@ impl<'a, S: GenericSim> RenderSteps<'a, S> {
         } else if false {
             me.render_e_z_field();
             me.render_b_xy_field();
-        } else {
+        } else if false {
             me.render_b();
             me.render_current();
+        } else {
+            me.render_m();
         }
         me.render_measurements();
         me.im
@@ -156,6 +158,11 @@ impl<'a, S: GenericSim> RenderSteps<'a, S> {
         self.render_vector_field(Rgb([0xff, 0xff, 0xff]), 1.0e-9, |cell| cell.b().xy());
     }
 
+    fn render_m(&mut self) {
+        self.render_scalar_field(1.0e5, false, 1, |cell| cell.mat().m().mag());
+        self.render_vector_field(Rgb([0xff, 0xff, 0xff]), 1.0e5, |cell| cell.mat().m().xy());
+    }
+
     fn render_vector_field<F: Fn(&Cell<mat::Static>) -> Vec2>(&mut self, color: Rgb<u8>, typical: Flt, measure: F) {
         let w = self.im.width();
         let h = self.im.height();
@@ -255,6 +262,12 @@ impl ImageRenderExt for RgbImage {
             pixelops::interpolate(left, right, left_weight*alpha)
         };
         drawing::draw_antialiased_line_segment_mut(self, i_start, i_end, color, interpolate_with_alpha);
+        if i_start != i_end
+            && (0..self.width() as i32).contains(&i_end.0)
+            && (0..self.height() as i32).contains(&i_end.1)
+        {
+            self.put_pixel(i_end.0 as _, i_end.1 as _, Rgb([0xff, 0, 0]));
+        }
         //drawing::draw_line_segment_mut(self, i_start, i_end, color);
     }
 }