Snapshot and decimate before rendering; add an explicit boundary layer

examples/ferromagnet.rs

@@ -1,14 +1,16 @@
 use coremem::{Driver, mat};
 
 fn main() {
-    let width = 201;
+    let width = 1300;
-    let height = 101;
+    let height = 1000;
-    let mut driver = Driver::new(width, height, 1e-3 /* feature size */)
+    let mut driver = Driver::new(width, height, 1e-5 /* feature size */)
         .with_y4m_renderer("ferromagnet.y4m")
+        .with_term_renderer()
+        .with_steps_per_frame(40)
         ;
 
     for y in 0..height {
-        for x in 50..60 {
+        for x in 500..600 {
             *driver.state.get_mut(x, y).mat_mut() = mat::Static::conductor(1.0e1).into();
         }
         // for x in 30..40 {
@@ -19,12 +21,12 @@ fn main() {
         //         *state.get_mut(x, y).mat_mut() = mat::Conductor { conductivity: 1.0e8 }.into();
         //     }
         // }
-        for x in 72..80 {
+        for x in 720..800 {
             *driver.state.get_mut(x, y).mat_mut() = mat::Static::conductor(1.0e1).into();
         }
     }
-    for y in 40..60 {
+    for y in 400..600 {
-        for x in 62..70 {
+        for x in 620..700 {
             *driver.state.get_mut(x, y).mat_mut() = mat::PiecewiseLinearFerromagnet::from_bh(&[
                 (  35.0, 0.0),
                 (  50.0, 0.250),
@@ -44,8 +46,8 @@ fn main() {
             //}.into();
         }
     }
+    driver.add_boundary(100);
 
-    //let mut renderer = render::NullRenderer;
     loop {
         //let imp = match state.step_no() {
         //    20..=60 => 1e6,
@@ -61,8 +63,8 @@ fn main() {
         // state.impulse_ey(50, 50, imp);
         // state.impulse_bz(20, 20, (imp / 3.0e8) as _);
         // state.impulse_bz(80, 20, (imp / 3.0e8) as _);
-        for y in 10..height-10 {
+        for y in 100..height-100 {
-            for x in 52..58 {
+            for x in 520..580 {
                 driver.state.impulse_ey(x, y, imp as _);
             }
         }

src/driver.rs

@@ -1,5 +1,6 @@
-use crate::sim::SimState;
+use crate::mat;
 use crate::render::{self, MultiRenderer, Renderer};
+use crate::sim::SimState;
 
 use std::path::PathBuf;
 use std::time::{Duration, SystemTime};
@@ -39,6 +40,26 @@ impl Driver {
         self.with_renderer(render::ColorTermRenderer)
     }
 
+    pub fn add_boundary(&mut self, thickness: u32) {
+        for inset in 0..thickness {
+            let depth = thickness - inset;
+            // TODO: tune a scalar multiplier on this value
+            let conductivity = (depth*depth) as f64;
+            for x in inset as usize..self.state.width() - inset as usize {
+                // left
+                *self.state.get_mut(x, inset as _).mat_mut() = mat::Static::conductor(conductivity).into();
+                // right
+                *self.state.get_mut(x, self.state.height() - 1 - inset as usize).mat_mut() = mat::Static::conductor(conductivity).into();
+            }
+            for y in inset as usize..self.state.height() - inset as usize {
+                // top
+                *self.state.get_mut(inset as _, y).mat_mut() = mat::Static::conductor(conductivity).into();
+                // bottom
+                *self.state.get_mut(self.state.width() - 1 - inset as usize, y).mat_mut() = mat::Static::conductor(conductivity).into();
+            }
+        }
+    }
+
     pub fn step(&mut self) {
         let start_time = SystemTime::now();
         if self.state.step_no() % self.steps_per_frame == 0 {

src/render.rs

@@ -1,6 +1,6 @@
 use ansi_term::Color::RGB;
 use crate::geom::Point;
-use crate::{Material as _, SimState};
+use crate::{Material as _, SimSnapshot, SimState};
 use decorum::R64;
 use image::{RgbImage, Rgb};
 use imageproc::{pixelops, drawing};
@@ -49,12 +49,12 @@ fn scale_vector(x: Point, typical_mag: f64) -> Point {
     x.with_mag(new_mag)
 }
 
-trait SimStateRenderExt {
+trait SimSnapshotRenderExt {
     fn to_image(&self) -> RgbImage;
     fn e_vector(&self, xidx: u32, yidx: u32, size: u32) -> Point;
 }
 
-impl SimStateRenderExt for SimState {
+impl SimSnapshotRenderExt for SimSnapshot {
     fn to_image(&self) -> RgbImage {
         let w = self.width().try_into().unwrap();
         let h = self.height().try_into().unwrap();
@@ -64,7 +64,7 @@ impl SimStateRenderExt for SimState {
             for x in 0..w {
                 let cell = self.get(x as usize, y as usize);
                 let r = scale_signed_to_u8(cell.mat().mz(), 100.0);
-                let b = scale_unsigned_to_u8(cell.mat().conductivity(), 3.0);
+                let b = scale_unsigned_to_u8(cell.mat().conductivity(), 10.0);
                 let g = scale_signed_to_u8(cell.bz(), 1.0e-4);
                 image.put_pixel(x, y, Rgb([r, g, b]));
             }
@@ -128,10 +128,10 @@ impl ImageRenderExt for RgbImage {
 }
 
 pub trait Renderer {
-    fn render(&mut self, state: &SimState) {
+    fn render(&mut self, state: &SimSnapshot) {
         self.render_with_image(state, &state.to_image());
     }
-    fn render_with_image(&mut self, state: &SimState, _im: &RgbImage) {
+    fn render_with_image(&mut self, state: &SimSnapshot, _im: &RgbImage) {
         self.render(state);
     }
 }
@@ -139,7 +139,7 @@ pub trait Renderer {
 pub struct NumericTermRenderer;
 
 impl Renderer for NumericTermRenderer {
-    fn render(&mut self, state: &SimState) {
+    fn render(&mut self, state: &SimSnapshot) {
         for y in 0..state.height() {
             for x in 0..state.width() {
                 let cell = state.get(x, y);
@@ -159,7 +159,7 @@ impl Renderer for NumericTermRenderer {
 pub struct ColorTermRenderer;
 
 impl Renderer for ColorTermRenderer {
-    fn render_with_image(&mut self, state: &SimState, im: &RgbImage) {
+    fn render_with_image(&mut self, state: &SimSnapshot, im: &RgbImage) {
         let square = "█";
         let buf: String = im
             .enumerate_rows()
@@ -189,7 +189,7 @@ impl Y4MRenderer {
 }
 
 impl Renderer for Y4MRenderer {
-    fn render_with_image(&mut self, _state: &SimState, im: &RgbImage) {
+    fn render_with_image(&mut self, _state: &SimSnapshot, im: &RgbImage) {
         if self.encoder.is_none() {
             let writer = File::create(&self.out_path).unwrap();
             self.encoder = Some(y4m::encode(im.width() as usize, im.height() as usize, y4m::Ratio::new(30, 1))
@@ -236,16 +236,24 @@ impl MultiRenderer {
         self.push(r);
         self
     }
+
+    pub fn render(&mut self, state: &SimState) {
+        //let max_width = 1980; //< TODO: make configurable
+        let max_width = 200;
+        let dec = (state.width() as u32 + max_width - 1) / max_width;
+        let snap = state.snapshot(dec);
+        Renderer::render(self, &snap);
+    }
 }
 
 impl Renderer for MultiRenderer {
-    fn render(&mut self, state: &SimState) {
+    fn render(&mut self, state: &SimSnapshot) {
         if self.renderers.len() != 0 {
             self.render_with_image(state, &state.to_image());
         }
     }
 
-    fn render_with_image(&mut self, state: &SimState, im: &RgbImage) {
+    fn render_with_image(&mut self, state: &SimSnapshot, im: &RgbImage) {
         for r in &mut self.renderers {
             r.render_with_image(state, im);
         }

src/sim.rs

@@ -1,11 +1,13 @@
 use crate::{consts};
 use crate::geom::Point;
-use crate::mat::{GenericMaterial, Material};
+use crate::mat::{self, GenericMaterial, Material};
 
 use decorum::R64;
 use ndarray::Array2;
 use rayon::prelude::*;
 
+pub type SimSnapshot = SimState<mat::Static>;
+
 #[derive(Default)]
 pub struct SimState<M=GenericMaterial> {
     cells: Array2<Cell<M>>,
@@ -70,6 +72,49 @@ impl<M: Material + Clone + Default + Send + Sync> SimState<M> {
 
         self.step_no += 1;
     }
+
+    /// dec = decimation, e.g. set to 3 to sample 3 x 3 grids
+    pub fn snapshot(&self, dec: u32) -> SimSnapshot {
+        let rows = self.height() as u32 / dec;
+        let cols = self.width() as u32 / dec;
+        let sample_area = (dec * dec) as f64;
+
+        let immute_self = &self;
+        let new_cells_vec: Vec<Vec<_>> = (0..rows).into_par_iter()
+            .map(|new_y| {
+                (0..cols).into_iter().map(move |new_x| {
+                    // sample dec x dec cells
+                    let mut cell: Cell<mat::Static> = Default::default();
+                    for y in new_y*dec..new_y*dec + dec {
+                        for x in new_x*dec..new_x*dec + dec {
+                            let sample = immute_self.get(x as _, y as _);
+                            cell.state.ex += sample.ex();
+                            cell.state.ey += sample.ey();
+                            cell.state.hz += sample.hz();
+                            cell.mat.conductivity += sample.mat.conductivity();
+                            cell.mat.mz += sample.mat.mz();
+                        }
+                    }
+                    cell.state.ex /= sample_area;
+                    cell.state.ey /= sample_area;
+                    cell.state.hz /= sample_area;
+                    cell.mat.conductivity /= sample_area;
+                    cell.mat.mz /= sample_area;
+                    cell
+                }).collect()
+            }).collect();
+
+        let mut new_cells: Array2<Cell<mat::Static>> = Array2::default((rows as _, cols as _));
+        for (row, mut new_row) in new_cells_vec.into_iter().enumerate() {
+            new_cells.row_mut(row).as_slice_mut().unwrap().swap_with_slice(&mut new_row[..]);
+        }
+
+        SimState {
+            cells: new_cells,
+            feature_size: self.feature_size,
+            step_no: self.step_no,
+        }
+    }
 }
 
 impl<M: Material> SimState<M> {
@@ -87,6 +132,10 @@ impl<M> SimState<M> {
         self.step_no
     }
 
+    pub fn feature_size(&self) -> f64 {
+        self.feature_size.into()
+    }
+
     pub fn impulse_ex(&mut self, x: usize, y: usize, ex: f64) {
         self.cells[[y, x]].state.ex += ex;
     }
@@ -227,13 +276,13 @@ impl<M: Material> Cell<M> {
         let sigma: R64 = self.mat.conductivity().into();
 
         // XXX not obvious that bz_to_hz is sensible
-        let delta_hz_y = self.hz() - self.bz_to_hz(up.bz());
+        // let delta_hz_y = self.hz() - self.bz_to_hz(up.bz());
-        // let delta_hz_y = self.hz() - up.hz();
+        let delta_hz_y = self.hz() - up.hz();
         let ex_rhs = self.state.ex*(ONE() - sigma/EPS0()*delta_t) + TWO()*delta_t/EPS0()*delta_hz_y/feature_size;
         let ex_next = ex_rhs / (ONE() + sigma/EPS0()*delta_t);
 
-        let delta_hz_x = self.hz() - self.bz_to_hz(left.bz());
+        // let delta_hz_x = self.hz() - self.bz_to_hz(left.bz());
-        // let delta_hz_x = self.hz() - left.hz();
+        let delta_hz_x = self.hz() - left.hz();
         let ey_rhs = self.state.ey*(ONE() - sigma/EPS0()*delta_t) - TWO()*delta_t/EPS0()*delta_hz_x/feature_size;
         let ey_next = ey_rhs / (ONE() + sigma/EPS0()*delta_t);