UPML boundary now applies to z dimension

examples/toroid25d.rs

@@ -7,23 +7,28 @@ fn main() {
     let from_m = |m| (m/feat_size) as u32;
     let m_to_um = |px| (px * 1e6) as u32;
     let to_m = |px| px as Flt * feat_size;
-    let width = 2500e-6;
+    let width = 2000e-6;
+    let depth = 500e-6;
     let conductor_inner_rad = 0e-6;
     let conductor_outer_rad = 200e-6;
     let ferro_inner_rad = 220e-6;
     let ferro_outer_rad = 300e-6;
-    let buffer = 200e-6;
+    let buffer = 100e-6;
     let peak_current = 100.0;
     let conductivity = 1.0e5;
     let half_width = width * 0.5;
+    let half_depth = depth * 0.5;
 
     let width_px = from_m(width);
-    let mut driver = Driver::new(width_px, width_px, feat_size);
+    let depth_px = from_m(depth);
+    let size_px = (width_px, width_px, depth_px).into();
+    let mut driver = Driver::new(size_px, feat_size);
     //driver.set_steps_per_frame(8);
     //driver.set_steps_per_frame(40);
     //driver.set_steps_per_frame(200);
-    driver.add_y4m_renderer(&*format!("toroid25d-flt{}-feat{}um-{}A--radii{}um-{}um-{}um.y4m",
+    driver.add_y4m_renderer(&*format!("toroid25d-flt{}-{}-feat{}um-{}A--radii{}um-{}um-{}um.y4m",
         std::mem::size_of::<Flt>() * 8,
+        size_px,
         m_to_um(feat_size),
         peak_current as u32,
         m_to_um(conductor_outer_rad),
@@ -37,50 +42,55 @@ fn main() {
         conductor_outer_rad)
     ));
     driver.add_measurement(meas::Magnetization(
-        (half_width + ferro_inner_rad + 2.0*feat_size, half_width, 0.0).into()
+        (half_width + ferro_inner_rad + 2.0*feat_size, half_width, half_depth).into()
     ));
     driver.add_measurement(meas::MagneticFlux(
-        (half_width + ferro_inner_rad + 2.0*feat_size, half_width, 0.0).into()
+        (half_width + ferro_inner_rad + 2.0*feat_size, half_width, half_depth).into()
     ));
     driver.add_measurement(meas::MagneticStrength(
-        (half_width + ferro_inner_rad + 2.0*feat_size, half_width, 0.0).into()
+        (half_width + ferro_inner_rad + 2.0*feat_size, half_width, half_depth).into()
     ));
     driver.add_measurement(meas::Magnetization(
-        (half_width + ferro_inner_rad + 1.0*feat_size, half_width, 0.0).into()
+        (half_width + ferro_inner_rad + 1.0*feat_size, half_width, half_depth).into()
     ));
     driver.add_measurement(meas::MagneticFlux(
-        (half_width + ferro_inner_rad + 1.0*feat_size, half_width, 0.0).into()
+        (half_width + ferro_inner_rad + 1.0*feat_size, half_width, half_depth).into()
     ));
     driver.add_measurement(meas::MagneticStrength(
-        (half_width + ferro_inner_rad + 1.0*feat_size, half_width, 0.0).into()
+        (half_width + ferro_inner_rad + 1.0*feat_size, half_width, half_depth).into()
     ));
     driver.add_measurement(meas::Magnetization(
-        (half_width + 0.5 * (ferro_inner_rad + ferro_outer_rad), half_width, 0.0).into()
+        (half_width + 0.5 * (ferro_inner_rad + ferro_outer_rad), half_width, half_depth).into()
     ));
     driver.add_measurement(meas::MagneticFlux(
-        (half_width + 0.5 * (ferro_inner_rad + ferro_outer_rad), half_width, 0.0).into()
+        (half_width + 0.5 * (ferro_inner_rad + ferro_outer_rad), half_width, half_depth).into()
     ));
     driver.add_measurement(meas::MagneticStrength(
-        (half_width + 0.5 * (ferro_inner_rad + ferro_outer_rad), half_width, 0.0).into()
+        (half_width + 0.5 * (ferro_inner_rad + ferro_outer_rad), half_width, half_depth).into()
     ));
 
     let center = Vec2::new(half_width as _, half_width as _);
 
+    for z_px in 0..depth_px {
         for y_px in 0..width_px {
             for x_px in 0..width_px { 
-            let loc = Coord::new(x_px, y_px, 0);
+                let loc = Coord::new(x_px, y_px, z_px);
                 let d = Vec2::new(to_m(x_px), to_m(y_px)) - center;
                 let r = d.mag();
                 if (conductor_inner_rad as _..conductor_outer_rad as _).contains(&r) {
                     *driver.state.get_mut(loc).mat_mut() = mat::Static::conductor(conductivity).into();
                 } else if (ferro_inner_rad as _..ferro_outer_rad as _).contains(&r) {
+                    //if (half_depth_px-5..half_depth_px+5).contains(z_px) {
                         *driver.state.get_mut(loc).mat_mut() = mat::db::ferroxcube_3r1();
+                    //}
                 }
             }
         }
-    let boundary = half_width - ferro_outer_rad - buffer;
+    }
-    println!("boundary: {}um", m_to_um(boundary));
+    let boundary_xy = half_width - ferro_outer_rad - buffer;
-    driver.add_upml_boundary(from_m(boundary));
+    println!("boundary: {}um", m_to_um(boundary_xy));
+    let boundary = Coord::new(from_m(boundary_xy), from_m(boundary_xy), 20);
+    driver.add_upml_boundary(boundary);
 
     loop {
         // let drive_current = peak_current * match driver.state.step_no() {
@@ -88,7 +98,7 @@ fn main() {
         //     3000..=4000 => -1.0,
         //     _ => 0.0,
         // };
-        let drive_current = peak_current;
+        let drive_current = peak_current * 1e-18;
         // J = \sigma*E = [Am^-2]
         // I = \sigma*E*Area
         // E = I / \sigma / Area

src/driver.rs

@@ -21,9 +21,9 @@ pub struct Driver {
 
 impl Driver {
     // TODO: allow depth
-    pub fn new(width: u32, height: u32, feature_size: Flt) -> Self {
+    pub fn new(size: Coord, feature_size: Flt) -> Self {
         Driver {
-            state: SimState::new(Coord::new(width, height, 1), feature_size),
+            state: SimState::new(size, feature_size),
             renderer: Default::default(),
             steps_per_frame: 1,
             time_spent_stepping: Default::default(),
@@ -75,34 +75,44 @@ impl Driver {
     //     }
     // }
 
-    // TODO: z-coordinate!
+    pub fn add_upml_boundary(&mut self, thickness: Coord) {
-    pub fn add_upml_boundary(&mut self, thickness: u32) {
         // Based on explanation here (slide 63): https://empossible.net/wp-content/uploads/2020/01/Lecture-The-Perfectly-Matched-Layer.pdf
         let h = self.state.height();
         let w = self.state.width();
-        for y in 0..h {
+        let d = self.state.depth();
-            for x in 0..w {
+        for z in 0..d {
-                let depth_x = if x < thickness {
+            let depth_z = if z < thickness.z() {
-                    thickness - x
+                thickness.z() - z
-                } else if x >= w - thickness {
+            } else if z >= d - thickness.z() {
-                    1 + x - (w - thickness)
+                1 + z - (d - thickness.z())
             } else {
                 0
             };
-                let depth_y = if y < thickness {
+            for y in 0..h {
-                    thickness - y
+                let depth_y = if y < thickness.y() {
-                } else if y >= h - thickness {
+                    thickness.y() - y
-                    1 + y - (h - thickness)
+                } else if y >= h - thickness.y() {
+                    1 + y - (h - thickness.y())
+                } else {
+                    0
+                };
+                for x in 0..w {
+                    let depth_x = if x < thickness.x() {
+                        thickness.x() - x
+                    } else if x >= w - thickness.x() {
+                        1 + x - (w - thickness.x())
                     } else {
                         0
                     };
 
-                if depth_x > 0 || depth_y > 0 {
+                    if depth_x > 0 || depth_y > 0 || depth_z > 0 {
                         let scale = 0.5 * consts::EPS0 / self.state.timestep();
-                    let cond_x = scale * (depth_x as Flt/thickness as Flt).powf(3.0);
+                        let cond_x = scale * (depth_x as Flt/thickness.x() as Flt).powf(3.0);
-                    let cond_y = scale * (depth_y as Flt/thickness as Flt).powf(3.0);
+                        let cond_y = scale * (depth_y as Flt/thickness.y() as Flt).powf(3.0);
-                    let conductor = mat::Static::anisotropic_conductor(Vec3::new(cond_x, cond_y, 0.0));
+                        let cond_z = scale * (depth_z as Flt/thickness.z() as Flt).powf(3.0);
-                    *self.state.get_mut((x, y, 0).into()).mat_mut() = conductor.into();
+                        let conductor = mat::Static::anisotropic_conductor(Vec3::new(cond_x, cond_y, cond_z));
+                        *self.state.get_mut((x, y, z).into()).mat_mut() = conductor.into();
+                    }
                 }
             }
         }

src/render.rs

@@ -56,13 +56,15 @@ struct RenderSteps<'a> {
     im: RgbImage, 
     sim: &'a dyn GenericSim,
     meas: &'a [Box<dyn AbstractMeasurement>],
+    /// Simulation z coordinate to sample
+    z: u32,
 }
 
 impl<'a> RenderSteps<'a> {
-    fn render(state: &'a dyn GenericSim, measurements: &'a [Box<dyn AbstractMeasurement>]) -> RgbImage {
+    fn render(state: &'a dyn GenericSim, measurements: &'a [Box<dyn AbstractMeasurement>], z: u32) -> RgbImage {
         let width = 768;
         let height = width * state.height() / state.width();
-        let mut me = Self::new(state, measurements, width, height);
+        let mut me = Self::new(state, measurements, width, height, z);
         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 {
@@ -75,11 +77,12 @@ impl<'a> RenderSteps<'a> {
         me.render_measurements();
         me.im
     }
-    fn new(sim: &'a dyn GenericSim, meas: &'a [Box<dyn AbstractMeasurement>], width: u32, height: u32) -> Self {
+    fn new(sim: &'a dyn GenericSim, meas: &'a [Box<dyn AbstractMeasurement>], width: u32, height: u32, z: u32) -> Self {
         RenderSteps {
             im: RgbImage::new(width, height),
             sim,
-            meas
+            meas,
+            z
         }
     }
 
@@ -88,7 +91,8 @@ impl<'a> RenderSteps<'a> {
         let x_m = x_prop * (self.sim.width() as Flt * self.sim.feature_size());
         let y_prop = y_px as Flt / self.im.height() as Flt;
         let y_m = y_prop * (self.sim.height() as Flt * self.sim.feature_size());
-        self.sim.sample(Vec3::new(x_m, y_m, 0.0))
+        let z_m = self.z as Flt * self.sim.feature_size();
+        self.sim.sample(Vec3::new(x_m, y_m, z_m))
     }
 
     //////////////   Ex/Ey/Bz configuration ////////////
@@ -212,7 +216,7 @@ impl ImageRenderExt for RgbImage {
 
 pub trait Renderer {
     fn render(&mut self, state: &dyn GenericSim, measurements: &[Box<dyn AbstractMeasurement>]) {
-        self.render_with_image(state, &RenderSteps::render(state, measurements));
+        self.render_with_image(state, &RenderSteps::render(state, measurements, state.depth() / 2));
     }
     fn render_with_image(&mut self, state: &dyn GenericSim, _im: &RgbImage) {
         self.render(state, &[]);
@@ -327,7 +331,7 @@ impl MultiRenderer {
 impl Renderer for MultiRenderer {
     fn render(&mut self, state: &dyn GenericSim, measurements: &[Box<dyn AbstractMeasurement>]) {
         if self.renderers.len() != 0 {
-            self.render_with_image(state, &RenderSteps::render(state, measurements));
+            self.render_with_image(state, &RenderSteps::render(state, measurements, state.depth() / 2));
         }
     }