From 1d303df409fe3eea0cfee259583b8fe7fa288424 Mon Sep 17 00:00:00 2001 From: Colin Date: Fri, 18 Sep 2020 22:22:54 -0700 Subject: [PATCH] Fix rendering to render (Bx,By,Ez). Also fix an infinite loop due to rounding errors in MHCurve --- src/driver.rs | 20 ++++++++++---- src/geom.rs | 7 +++++ src/mat.rs | 37 ++++++++++++++++++++++--- src/render.rs | 75 ++++++++++++++++++++++++++++++++++----------------- src/sim.rs | 5 ++++ src/vec3.rs | 6 ++++- 6 files changed, 116 insertions(+), 34 deletions(-) diff --git a/src/driver.rs b/src/driver.rs index 9a5d1db..d62f562 100644 --- a/src/driver.rs +++ b/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>, + 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); } } } diff --git a/src/geom.rs b/src/geom.rs index ced4848..0bdd9d5 100644 --- a/src/geom.rs +++ b/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() diff --git a/src/mat.rs b/src/mat.rs index 41abff7..5c80e95 100644 --- a/src/mat.rs +++ b/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); + } } diff --git a/src/render.rs b/src/render.rs index b5c1cd7..08fb580 100644 --- a/src/render.rs +++ b/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]) -> 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) -> Point>(&mut self, color: Rgb, 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) -> 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 } } diff --git a/src/sim.rs b/src/sim.rs index 68f89c2..0ac6766 100644 --- a/src/sim.rs +++ b/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 SimState { 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; diff --git a/src/vec3.rs b/src/vec3.rs index f0dcc9c..3a88f0f 100644 --- a/src/vec3.rs +++ b/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 for Vec3 { impl DivAssign for Vec3 { fn div_assign(&mut self, other: Flt) { - *self *= (1.0 / other); + *self *= 1.0 / other; } }