Allow conductivity to be specified as a vector

This is necessary to implement UPML boundary conditions

examples/toroid25d.rs

@@ -3,20 +3,29 @@ use coremem::geom::Point;
 
 fn main() {
     coremem::init_logging();
-    let width = 700;
+    let feat_size = 20e-6; // feature size
-    let buffer = 80;
+    let from_m = |m| (m/feat_size) as u32;
-    let feat_size = 1e-5; // feature size
+    let to_um = |px| px * (feat_size*1e6) as u32;
+    let width = from_m(8000e-6);
+    let conductor_inner_rad = from_m(0e-6);
+    let conductor_outer_rad = from_m(200e-6);
+    let ferro_inner_rad = from_m(200e-6);
+    let ferro_outer_rad = from_m(300e-6);
+    let buffer = from_m(800e-6);
     let peak_current = 2.5e6;
     let conductivity = 1.0e5;
-    let ferro_inner_rad = 60;
-    let ferro_outer_rad = 150;
-    let conductor_inner_rad = 0;
-    let conductor_outer_rad = 50;
     let half_width = width / 2;
     let mut driver = Driver::new(width, width, feat_size);
-    driver.set_steps_per_frame(8);
+    //driver.set_steps_per_frame(8);
     //driver.set_steps_per_frame(40);
-    driver.add_y4m_renderer("toroid25d.y4m");
+    driver.set_steps_per_frame(200);
+    driver.add_y4m_renderer(&*format!("-toroid25d-feat{}um-{}A--radii{}um-{}um-{}um.y4m",
+        (feat_size * 1e6) as u32,
+        peak_current as u32,
+        to_um(conductor_outer_rad),
+        to_um(ferro_inner_rad),
+        to_um(ferro_outer_rad),
+    ));
     // driver.add_term_renderer();
     driver.add_measurement(meas::Label(format!("Conductivity: {}, Imax: {:.2e}", conductivity, peak_current)));
     driver.add_measurement(meas::Current(half_width + (conductor_inner_rad + conductor_outer_rad) / 2, half_width));

src/driver.rs

@@ -12,6 +12,7 @@ pub struct Driver {
     renderer: MultiRenderer,
     steps_per_frame: u64,
     time_spent_stepping: Duration,
+    time_spent_rendering: Duration,
     measurements: Vec<Box<dyn AbstractMeasurement>>,
     start_time: SystemTime,
 }
@@ -23,6 +24,7 @@ impl Driver {
             renderer: Default::default(),
             steps_per_frame: 1,
             time_spent_stepping: Default::default(),
+            time_spent_rendering: Default::default(),
             measurements: vec![Box::new(meas::Time), Box::new(meas::Meta), Box::new(meas::Energy)],
             start_time: SystemTime::now(),
         }
@@ -71,22 +73,34 @@ 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");
+            let start_time = SystemTime::now();
             self.renderer.render(&self.state, &*self.measurements);
+            self.time_spent_rendering += start_time.elapsed().unwrap();
             trace!("render end");
         }
         trace!("step begin");
+        let start_time = SystemTime::now();
         self.state.step();
-        trace!("step end");
         self.time_spent_stepping += start_time.elapsed().unwrap();
+        trace!("step end");
         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 render_time = self.time_spent_rendering.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);
+            let sim_time = self.state.time();
+            info!(
+                "t={:.2e} frame {:06} fps: {:6.2} (sim: {:.1}s, render: {:.1}s, other: {:.1}s)",
+                sim_time,
+                step,
+                fps,
+                driver_time,
+                render_time,
+                overall_time - driver_time - render_time
+            );
         }
     }
 }

src/mat.rs

@@ -10,8 +10,8 @@ use std::cmp::Ordering;
 pub trait Material {
     /// Return \sigma, the electrical conductivity.
     /// For a vacuum, this is zero. For a perfect conductor, \inf.
-    fn conductivity(&self) -> Flt {
+    fn conductivity(&self) -> Vec3 {
-        0.0
+        Default::default()
     }
     /// Return the magnetization.
     fn m(&self) -> Vec3 {
@@ -26,21 +26,21 @@ pub trait Material {
 
 #[derive(Clone, Default)]
 pub struct Static {
-    pub conductivity: Flt,
+    pub conductivity: Vec3,
     pub m: Vec3,
 }
 
 impl Static {
     pub fn conductor(conductivity: Flt) -> Self {
         Self {
-            conductivity,
+            conductivity: Vec3::unit() * conductivity,
             ..Default::default()
         }
     }
 }
 
 impl Material for Static {
-    fn conductivity(&self) -> Flt {
+    fn conductivity(&self) -> Vec3 {
         self.conductivity
     }
     fn m(&self) -> Vec3 {
@@ -140,7 +140,7 @@ impl MHCurve {
 pub struct PiecewiseLinearFerromagnet {
     /// Instantaneous magnetization in the z direction.
     pub m: Vec3,
-    pub conductivity: Flt,
+    pub conductivity: Vec3,
     mh_curve_x: MHCurve,
     mh_curve_y: MHCurve,
     mh_curve_z: MHCurve,
@@ -157,7 +157,7 @@ impl PiecewiseLinearFerromagnet {
 
         Self {
             m: Vec3::default(),
-            conductivity: 0.0,
+            conductivity: Default::default(),
             mh_curve_x: MHCurve::new(&*mh_points),
             mh_curve_y: MHCurve::new(&*mh_points),
             mh_curve_z: MHCurve::new(&*mh_points),
@@ -169,7 +169,7 @@ impl Material for PiecewiseLinearFerromagnet {
     fn m(&self) -> Vec3 {
         self.m
     }
-    fn conductivity(&self) -> Flt {
+    fn conductivity(&self) -> Vec3 {
         self.conductivity
     }
     fn step_h(&mut self, context: &CellState, delta_b: Vec3) -> Vec3 {
@@ -223,7 +223,7 @@ pub mod db {
             (  50.0, 0.340),
             (   0.0, 0.325),
         ]);
-        m.conductivity = 1e3;
+        m.conductivity = Vec3::unit() * 1e3; // XXX colin: should be 1e-3!!?
         m
     }
 }

src/render.rs

@@ -61,7 +61,7 @@ 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_scalar_field(10.0,  false, 2, |cell| cell.mat().conductivity());
+        me.render_scalar_field(10.0,  false, 2, |cell| cell.mat().conductivity().mag());
         me.render_scalar_field(100.0, true,  0, |cell| cell.mat().m().mag());
         if false {
             me.render_b_z_field();
@@ -91,7 +91,7 @@ impl<'a> RenderSteps<'a> {
         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()
+            cell.e().elem_mul(cell.mat().conductivity()).xy()
         });
     }
 

src/sim.rs

@@ -219,7 +219,7 @@ impl<M: Material> Cell<M> {
 
     pub fn current(&self) -> Vec3 {
         let conductivity = self.mat.conductivity();
-        self.e()*conductivity
+        self.e().elem_mul(conductivity)
     }
 
     fn impulse_bz(&mut self, delta_bz: Flt) {
@@ -288,7 +288,7 @@ impl<M: Material> Cell<M> {
         
         use consts::real::{EPS0, ONE, TWO};
 
-        let sigma: Real = self.mat.conductivity().into();
+        let sigma: Vec3 = self.mat.conductivity();
         let inv_feature_size = Real::from_inner(1.0) / feature_size;
 
         let delta_hz_y = inv_feature_size * (self.hz() - up.hz());
@@ -300,8 +300,8 @@ impl<M: Material> Cell<M> {
         let nabla_h = Vec3::from((delta_hz_y - delta_hy_z, delta_hx_z - delta_hz_x, delta_hy_x - delta_hx_y));
         let delta_t_eps = delta_t/EPS0();
 
-        let rhs = self.state.e() * (ONE() - sigma*delta_t_eps) + nabla_h * (TWO()*delta_t_eps);
+        let rhs = self.state.e().elem_mul(Vec3::unit() - sigma*delta_t_eps) + nabla_h * (TWO()*delta_t_eps);
-        let e_next = rhs / (ONE() + sigma*delta_t_eps);
+        let e_next = rhs.elem_div(Vec3::unit() + sigma*delta_t_eps);
 
 
         Cell {

src/vec3.rs

@@ -40,6 +40,24 @@ impl Vec3 {
         let Vec3 { x, y, z } = *self;
         (x*x + y*y + z*z).into_inner().sqrt()
     }
+
+    /// Perform element-wise multiplication with `other`.
+    pub fn elem_mul(&self, other: Self) -> Self {
+        Self {
+            x: self.x * other.x,
+            y: self.y * other.y,
+            z: self.z * other.z,
+        }
+    }
+
+    /// Perform element-wise division with `other`.
+    pub fn elem_div(&self, other: Self) -> Self {
+        Self {
+            x: self.x / other.x,
+            y: self.y / other.y,
+            z: self.z / other.z,
+        }
+    }
 }
 
 impl From<(Real, Real, Real)> for Vec3 {