Fix rendering to render (Bx,By,Ez). Also fix an infinite loop due to rounding errors in MHCurve

src/driver.rs

@@ -3,26 +3,28 @@ use crate::meas::{self, AbstractMeasurement};
 use crate::render::{self, MultiRenderer, Renderer};
 use crate::sim::SimState;
 
-use log::info;
+use log::{info, trace};
 use std::path::PathBuf;
 use std::time::{Duration, SystemTime};
 
-#[derive(Default)]
 pub struct Driver {
     pub state: SimState,
     renderer: MultiRenderer,
     steps_per_frame: u64,
     time_spent_stepping: Duration,
     measurements: Vec<Box<dyn AbstractMeasurement>>,
+    start_time: SystemTime,
 }
 
 impl Driver {
     pub fn new(width: u32, height: u32, feature_size: Flt) -> Self {
         Driver {
             state: SimState::new(width, height, feature_size),
+            renderer: Default::default(),
             steps_per_frame: 1,
+            time_spent_stepping: Default::default(),
             measurements: vec![Box::new(meas::Time), Box::new(meas::Meta), Box::new(meas::Energy)],
-            ..Default::default()
+            start_time: SystemTime::now(),
         }
     }
 
@@ -71,12 +73,20 @@ impl Driver {
     pub fn step(&mut self) {
         let start_time = SystemTime::now();
         if self.state.step_no() % self.steps_per_frame == 0 {
+            trace!("render begin");
             self.renderer.render(&self.state, &*self.measurements);
+            trace!("render end");
         }
+        trace!("step begin");
         self.state.step();
+        trace!("step end");
         self.time_spent_stepping += start_time.elapsed().unwrap();
-        if self.state.step_no() % (10*self.steps_per_frame) == 0 {
-            info!("fps: {:6.2}", (self.state.step_no() as Flt) / (self.time_spent_stepping.as_secs_f64() as Flt));
+        let step = self.state.step_no();
+        if step % (10*self.steps_per_frame) == 0 {
+            let driver_time = self.time_spent_stepping.as_secs_f64() as Flt;
+            let overall_time = self.start_time.elapsed().unwrap().as_secs_f64() as Flt;
+            let fps = (self.state.step_no() as Flt) / overall_time;
+            info!("frame {:6}/fps: {:6.2} (driver: {:.1}s, rest: {:.1}s)", step, fps, driver_time, overall_time - driver_time);
         }
     }
 }

src/geom.rs

@@ -140,6 +140,13 @@ impl Line {
         Point::new(x, self.y(x))
     }
 
+    pub fn from(&self) -> Point {
+        self.from
+    }
+    pub fn to(&self) -> Point {
+        self.to
+    }
+
     pub fn x_intercept(&self) -> Flt {
         let t = (-self.from.y) / (self.to.y - self.from.y);
         (self.from.x + t * (self.to.x - self.from.x)).into()

src/mat.rs

@@ -2,6 +2,7 @@ use crate::{CellState, consts};
 use crate::flt::Flt;
 use crate::geom::{Line, Point, Polygon};
 use crate::vec3::Vec3;
+use log::{debug, trace};
 use enum_dispatch::enum_dispatch;
 use std::cmp::Ordering;
 
@@ -93,10 +94,20 @@ impl MHCurve {
                 }
             }
         };
+        trace!("active segment: {:?}", active_segment);
 
         // Find some m(h) on the active_segment such that sum(h) = h + m(h) = target_hm
         let sum_h = active_segment + Line::new(Point::new(0.0, 0.0), Point::new(1.0, 1.0));
-        let new_h = sum_h.move_toward_y_unclamped(h, target_hm);
+        trace!("sum_h: {:?}", sum_h);
+        let new_h = if sum_h.to().y() != sum_h.from().y() {
+            sum_h.move_toward_y_unclamped(h, target_hm)
+        } else {
+            // avoid a division-by-zero.
+            // We could be anywhere along this line, but we prefer the endpoint
+            // so as to escape out of any permanent loops
+            active_segment.to().x()
+        };
+        trace!("new_h: {}", new_h);
 
         if sum_h.contains_x(new_h) {
             // the segment contains a point with the target H+M
@@ -107,7 +118,9 @@ impl MHCurve {
         }
     }
     fn move_to(&self, mut h: Flt, mut m: Flt, target_hm: Flt) -> (Flt, Flt) {
+        let mut i = 0;
         loop {
+            i += 1;
             match self.step_toward(h, m, target_hm) {
                 Ok((x, y)) => break (x, y),
                 Err((x, y)) => {
@@ -115,6 +128,10 @@ impl MHCurve {
                     m = y;
                 },
             }
+            if i % 256 == 0 {
+                debug!("unusually high iteration count without converging: {}. args: {}, {}, {}", i, h, m, target_hm);
+                panic!("die");
+            }
         }
     }
 }
@@ -156,6 +173,7 @@ impl Material for PiecewiseLinearFerromagnet {
         self.conductivity
     }
     fn step_h(&mut self, context: &CellState, delta_b: Vec3) -> Vec3 {
+        trace!("step_h enter");
         let (h, m) = (context.h(), self.m);
         let target_hm = h + m + delta_b/consts::real::MU0();
 
@@ -165,7 +183,9 @@ impl Material for PiecewiseLinearFerromagnet {
         let (hz, mz) = self.mh_curve_x.move_to(h.z(), m.z(), target_hm.z());
 
         self.m = Vec3::new(mx, my, mz);
-        Vec3::new(hx, hy, hz)
+        let ret = Vec3::new(hx, hy, hz);
+        trace!("step_h end");
+        ret
     }
 }
 
@@ -188,6 +208,9 @@ pub mod db {
     use super::*;
     /// https://www.ferroxcube.com/upload/media/product/file/MDS/3r1.pdf
     pub fn ferroxcube_3r1() -> GenericMaterial {
+        ferroxcube_3r1_concrete().into()
+    }
+    pub fn ferroxcube_3r1_concrete() -> PiecewiseLinearFerromagnet {
         let mut m = PiecewiseLinearFerromagnet::from_bh(&[
             (  35.0, 0.0),
             (  50.0, 0.250),
@@ -201,7 +224,7 @@ pub mod db {
             (   0.0, 0.325),
         ]);
         m.conductivity = 1e3;
-        m.into()
+        m
     }
 }
 
@@ -316,4 +339,12 @@ mod test {
         // Falling from interior to middle
         assert_move_to_symmetric(28.0, 130.0, 130.0, (5.0, 125.0));
     }
+
+    /// Float rounding would cause `inf`s, which manifested as infinite looping.
+    #[test]
+    fn regression_no_convergence_3r1() {
+        let mat = db::ferroxcube_3r1_concrete();
+        let curve = mat.mh_curve_x;
+        curve.move_to(-202.04596, -278400.53, -278748.66);
+    }
 }

src/render.rs

@@ -5,6 +5,7 @@ use crate::mat;
 use crate::sim::Cell;
 use crate::meas::AbstractMeasurement;
 use font8x8::{BASIC_FONTS, GREEK_FONTS, UnicodeFonts as _};
+use log::trace;
 use image::{RgbImage, Rgb};
 use imageproc::{pixelops, drawing};
 use std::fs::File;
@@ -60,8 +61,15 @@ struct RenderSteps<'a> {
 impl<'a> RenderSteps<'a> {
     fn render(state: &'a SimSnapshot, measurements: &'a [Box<dyn AbstractMeasurement>]) -> RgbImage {
         let mut me = Self::new(state, measurements);
-        me.render_b_field();
-        me.render_e_field();
+        me.render_scalar_field(10.0,  false, 2, |cell| cell.mat().conductivity());
+        me.render_scalar_field(100.0, true,  0, |cell| cell.mat().m().mag());
+        if false {
+            me.render_b_z_field();
+            me.render_e_xy_field();
+        } else {
+            me.render_e_z_field();
+            me.render_b_xy_field();
+        }
         me.render_measurements();
         me.im
     }
@@ -75,19 +83,26 @@ impl<'a> RenderSteps<'a> {
         }
     }
 
-    fn render_b_field(&mut self) {
-        let w = self.sim.width();
-        let h = self.sim.height();
-        for y in 0..h {
-            for x in 0..w {
-                let cell = self.sim.get(x, y);
-                let r = scale_signed_to_u8(cell.mat().m().z(), 100.0);
-                let b = scale_unsigned_to_u8(cell.mat().conductivity(), 10.0);
-                let g = scale_signed_to_u8(cell.b().z(), 1.0e-4);
-                self.im.put_pixel(x, y, Rgb([r, g, b]));
-            }
-        }
+    //////////////   Ex/Ey/Bz configuration ////////////
+    fn render_b_z_field(&mut self) {
+        self.render_scalar_field(1.0e-4, true,  1, |cell| cell.b().z());
     }
+    fn render_e_xy_field(&mut self) {
+        self.render_vector_field(Rgb([0xff, 0xff, 0xff]), 100.0, |cell| cell.e().xy());
+        // current
+        self.render_vector_field(Rgb([0x00, 0xa0, 0x30]), 1.0e-12, |cell| {
+            cell.e().xy()*cell.mat().conductivity()
+        });
+    }
+
+    //////////////   Bx/By/Ez configuration ////////////
+    fn render_e_z_field(&mut self) {
+        self.render_scalar_field(100.0, true,  1, |cell| cell.e().z());
+    }
+    fn render_b_xy_field(&mut self) {
+        self.render_vector_field(Rgb([0xff, 0xff, 0xff]), 1.0e-4, |cell| cell.b().xy());
+    }
+
     fn render_vector_field<F: Fn(&Cell<mat::Static>) -> Point>(&mut self, color: Rgb<u8>, typical: Flt, measure: F) {
         let w = self.sim.width();
         let h = self.sim.height();
@@ -105,16 +120,23 @@ impl<'a> RenderSteps<'a> {
             }
         }
     }
-    fn render_e_field(&mut self) {
-        self.render_vector_field(Rgb([0xff, 0xff, 0xff]), 100.0, |cell| cell.e().xy());
-        // current
-        self.render_vector_field(Rgb([0x00, 0xa0, 0x30]), 1.0e-12, |cell| {
-            //if cell.mat().conductivity() >= 1.0e3 {
-                cell.e().xy()*cell.mat().conductivity()
-            //} else {
-            //    Default::default()
-            //}
-        });
+    fn render_scalar_field<F: Fn(&Cell<mat::Static>) -> Flt>(&mut self, typical: Flt, signed: bool, slot: u32, measure: F) {
+        let w = self.sim.width();
+        let h = self.sim.height();
+        for y in 0..h {
+            for x in 0..w {
+                let cell = self.sim.get(x, y);
+                let value = measure(cell);
+                let scaled = if signed {
+                    scale_signed_to_u8(value, typical)
+                } else {
+                    scale_unsigned_to_u8(value, typical)
+                };
+                let mut p = *self.im.get_pixel(x, y);
+                p.0[slot as usize] = scaled;
+                self.im.put_pixel(x, y, p);
+            }
+        }
     }
     fn render_measurements(&mut self) {
         for (meas_no, m) in self.meas.iter().enumerate() {
@@ -269,7 +291,10 @@ impl Renderer for Y4MRenderer {
 
         let frame = y4m::Frame::new([&*pix_y, &*pix_u, &*pix_v], None);
         let enc = self.encoder.as_mut().unwrap();
-        enc.write_frame(&frame).unwrap()
+        trace!("write_frame begin");
+        let ret = enc.write_frame(&frame).unwrap();
+        trace!("write_frame end");
+        ret
     }
 }
 

src/sim.rs

@@ -2,6 +2,7 @@ use crate::{flt::{Flt, Real}, consts};
 use crate::geom::Point;
 use crate::mat::{self, GenericMaterial, Material};
 use crate::vec3::Vec3;
+use log::trace;
 
 use ndarray::{Array2, s, Zip};
 use ndarray::parallel::prelude::*;
@@ -39,20 +40,24 @@ impl<M: Material + Clone + Default + Send + Sync> SimState<M> {
         let mut scratch_cells = std::mem::replace(&mut self.scratch, Default::default());
 
         // first advance all the magnetic fields
+        trace!("step_h begin");
         Zip::from(self.cells.slice(s![0..h-1, 0..w-1]))
             .and(self.cells.slice(s![0..h-1, 1..w]))
             .and(self.cells.slice(s![1..h, 0..w-1]))
             .par_apply_assign_into(scratch_cells.slice_mut(s![0..h-1, 0..w-1]), |cell, right_cell, down_cell| {
                 cell.clone().step_h(right_cell, down_cell, half_time_step.into(), feature_size.into())
             });
+        trace!("step_h end");
 
         // now advance electic fields
+        trace!("step_e begin");
         Zip::from(scratch_cells.slice(s![1..h, 1..w]))
             .and(scratch_cells.slice(s![1..h, 0..w-1]))
             .and(scratch_cells.slice(s![0..h-1, 1..w]))
             .par_apply_assign_into(self.cells.slice_mut(s![1..h, 1..w]), |cell, left_cell, up_cell| {
                 cell.clone().step_e(left_cell, up_cell, half_time_step.into(), feature_size.into())
             });
+        trace!("step_e end");
 
         self.scratch = scratch_cells;
 

src/vec3.rs

@@ -36,6 +36,10 @@ impl Vec3 {
     pub fn xy(&self) -> Point {
         Point::new(self.x(), self.y())
     }
+    pub fn mag(&self) -> Flt {
+        let Vec3 { x, y, z } = *self;
+        (x*x + y*y + z*z).into_inner().sqrt()
+    }
 }
 
 impl From<(Real, Real, Real)> for Vec3 {
@@ -134,7 +138,7 @@ impl Div<Real> for Vec3 {
 
 impl DivAssign<Flt> for Vec3 {
     fn div_assign(&mut self, other: Flt) {
-        *self *= (1.0 / other);
+        *self *= 1.0 / other;
     }
 }