diff --git a/examples/toroid25d.rs b/examples/toroid25d.rs index d76fad5..b610e9d 100644 --- a/examples/toroid25d.rs +++ b/examples/toroid25d.rs @@ -1,106 +1,119 @@ use coremem::{Driver, Flt, mat, meas}; use coremem::geom::{Coord, CylinderZ, Index, Meters, Vec2, Vec3}; use coremem::stim::{Stimulus, Sinusoid}; +use log::trace; fn main() { coremem::init_logging(); - let feat_size = 10e-6; // feature size - 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 = 1000e-6; - let depth = 600e-6; - let conductor_inner_rad = 0e-6; - let conductor_outer_rad = 19e-6; - let ferro_inner_rad = 100e-6; - let ferro_outer_rad = 200e-6; - let buffer = 50e-6; - let peak_current = 2e3; - let current_duration = 1e-9; // half-wavelength of the sine wave - let conductivity = 1.0e5; - let half_width = width * 0.5; - let half_depth = depth * 0.5; + for p in 0..20 { + let ferro_depth = 10e-6 * (20-p) as Flt; + let feat_size = 10e-6; // feature size + 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 depth = 200e-6; + let conductor_inner_rad = 0e-6; + let conductor_outer_rad = 19e-6; + let ferro_inner_rad = 100e-6; + let ferro_outer_rad = 200e-6; + //let ferro_depth = 10e-6; + let buffer = 250e-6; + let peak_current = 2e3; + let current_duration = 1e-9; // half-wavelength of the sine wave + let conductivity = 1.0e5; + let half_width = width * 0.5; + let half_depth = depth * 0.5; - let width_px = from_m(width); - let depth_px = from_m(depth); - let size_px = Index((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-{:.1e}A-{:.1e}s--radii{}um-{}um-{}um.y4m", - std::mem::size_of::() * 8, - *size_px, - m_to_um(feat_size), - peak_current, - current_duration, - m_to_um(conductor_outer_rad), - m_to_um(ferro_inner_rad), - m_to_um(ferro_outer_rad), - )); - let conductor_region = CylinderZ::new( - Vec2::new(half_width, half_width), - conductor_outer_rad); - // driver.add_term_renderer(); - driver.add_measurement(meas::Label(format!("Conductivity: {}, Imax: {:.2e}", conductivity, peak_current))); - driver.add_measurement(meas::Current(conductor_region.clone())); - driver.add_measurement(meas::Magnetization( - Meters((half_width + ferro_inner_rad + 2.0*feat_size, half_width, half_depth).into()) - )); - driver.add_measurement(meas::MagneticFlux( - Meters((half_width + ferro_inner_rad + 2.0*feat_size, half_width, half_depth).into()) - )); - driver.add_measurement(meas::MagneticStrength( - Meters((half_width + ferro_inner_rad + 2.0*feat_size, half_width, half_depth).into()) - )); - driver.add_measurement(meas::Magnetization( - Meters((half_width + ferro_inner_rad + 1.0*feat_size, half_width, half_depth).into()) - )); - driver.add_measurement(meas::MagneticFlux( - Meters((half_width + ferro_inner_rad + 1.0*feat_size, half_width, half_depth).into()) - )); - driver.add_measurement(meas::MagneticStrength( - Meters((half_width + ferro_inner_rad + 1.0*feat_size, half_width, half_depth).into()) - )); - driver.add_measurement(meas::Magnetization( - Meters((half_width + 0.5 * (ferro_inner_rad + ferro_outer_rad), half_width, half_depth).into()) - )); - driver.add_measurement(meas::MagneticFlux( - Meters((half_width + 0.5 * (ferro_inner_rad + ferro_outer_rad), half_width, half_depth).into()) - )); - driver.add_measurement(meas::MagneticStrength( - Meters((half_width + 0.5 * (ferro_inner_rad + ferro_outer_rad), half_width, half_depth).into()) - )); + let duration = 1e-9; - let center = Vec2::new(half_width as _, half_width as _); + let width_px = from_m(width); + let depth_px = from_m(depth); + let size_px = Index((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(80); + driver.set_steps_per_frame(120); + //driver.set_steps_per_frame(200); + driver.add_y4m_renderer(&*format!("toroid25d.5-flt{}-{}-feat{}um-{:.1e}A-{:.1e}s--radii{}um-{}um-{}um-{}um.y4m", + std::mem::size_of::() * 8, + *size_px, + m_to_um(feat_size), + peak_current, + current_duration, + m_to_um(conductor_outer_rad), + m_to_um(ferro_inner_rad), + m_to_um(ferro_outer_rad), + m_to_um(ferro_depth), + )); + let conductor_region = CylinderZ::new( + Vec2::new(half_width, half_width), + conductor_outer_rad); + // driver.add_term_renderer(); + driver.add_measurement(meas::Label(format!("Conductivity: {}, Imax: {:.2e}", conductivity, peak_current))); + driver.add_measurement(meas::Current(conductor_region.clone())); + driver.add_measurement(meas::Magnetization( + Meters((half_width + ferro_inner_rad + 2.0*feat_size, half_width, half_depth).into()) + )); + driver.add_measurement(meas::MagneticFlux( + Meters((half_width + ferro_inner_rad + 2.0*feat_size, half_width, half_depth).into()) + )); + driver.add_measurement(meas::MagneticStrength( + Meters((half_width + ferro_inner_rad + 2.0*feat_size, half_width, half_depth).into()) + )); + driver.add_measurement(meas::Magnetization( + Meters((half_width + ferro_inner_rad + 1.0*feat_size, half_width, half_depth).into()) + )); + driver.add_measurement(meas::MagneticFlux( + Meters((half_width + ferro_inner_rad + 1.0*feat_size, half_width, half_depth).into()) + )); + driver.add_measurement(meas::MagneticStrength( + Meters((half_width + ferro_inner_rad + 1.0*feat_size, half_width, half_depth).into()) + )); + driver.add_measurement(meas::Magnetization( + Meters((half_width + 0.5 * (ferro_inner_rad + ferro_outer_rad), half_width, half_depth).into()) + )); + driver.add_measurement(meas::MagneticFlux( + Meters((half_width + 0.5 * (ferro_inner_rad + ferro_outer_rad), half_width, half_depth).into()) + )); + driver.add_measurement(meas::MagneticStrength( + Meters((half_width + 0.5 * (ferro_inner_rad + ferro_outer_rad), half_width, half_depth).into()) + )); - for z_px in 0..depth_px { - for y_px in 0..width_px { - for x_px in 0..width_px { - let loc = Index((x_px, y_px, z_px).into()); - 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 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 = Index((x_px, y_px, z_px).into()); + 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) { + let half_depth_px = from_m(half_depth); + let ferro_depth_px = from_m(ferro_depth); + if (half_depth_px-ferro_depth_px/2 .. half_depth_px+(ferro_depth_px+1)/2).contains(&z_px) { + trace!("placing ferro at {:?}", loc); + *driver.state.get_mut(loc).mat_mut() = mat::db::ferroxcube_3r1(); + } + } } } } - } - let boundary_xy = half_width - ferro_outer_rad - buffer; - println!("boundary: {}um", m_to_um(boundary_xy)); - let boundary = Index((from_m(boundary_xy), from_m(boundary_xy), 0).into()); - driver.add_upml_boundary(boundary); + let boundary_xy = half_width - ferro_outer_rad - buffer; + println!("boundary: {}um", m_to_um(boundary_xy)); + let boundary = Index((from_m(boundary_xy), from_m(boundary_xy), 0).into()); + driver.add_upml_boundary(boundary); - driver.add_stimulus(Stimulus::new( - conductor_region.clone(), - Sinusoid::from_wavelength(Vec3::new(0.0, 0.0, peak_current), current_duration * 2.0 -))); + driver.add_stimulus(Stimulus::new( + conductor_region.clone(), + Sinusoid::from_wavelength(Vec3::new(0.0, 0.0, peak_current), current_duration * 2.0 + ))); - loop { - driver.step(); + while driver.dyn_state().time() < duration { + driver.step(); + } } } diff --git a/src/driver.rs b/src/driver.rs index 41e54a9..c9a5441 100644 --- a/src/driver.rs +++ b/src/driver.rs @@ -137,7 +137,7 @@ impl Driver { 0.0 }; let cond = Vec3::new(cond_x, cond_y, cond_z); - trace!("cond: {:?}", cond); + // trace!("cond: {:?}", cond); let conductor = mat::Static::anisotropic_conductor(cond); *self.state.get_mut(Index(Vec3u::new(x, y, z))).mat_mut() = conductor.into(); } diff --git a/src/geom/units.rs b/src/geom/units.rs index fd31d46..f204f62 100644 --- a/src/geom/units.rs +++ b/src/geom/units.rs @@ -9,7 +9,7 @@ pub trait Coord: Copy + Clone { fn from_index(other: Index, feature_size: Flt) -> Self; } -#[derive(Copy, Clone)] +#[derive(Copy, Clone, Debug)] pub struct Meters(pub Vec3); impl Coord for Meters { @@ -34,7 +34,7 @@ impl Deref for Meters { } } -#[derive(Copy, Clone)] +#[derive(Copy, Clone, Debug)] pub struct Index(pub Vec3u); impl Coord for Index { diff --git a/src/geom/vec.rs b/src/geom/vec.rs index fab4e83..5508cac 100644 --- a/src/geom/vec.rs +++ b/src/geom/vec.rs @@ -199,7 +199,7 @@ impl From<(Real, Real, Real)> for Vec3 { impl From for Vec3 { fn from(v: Vec3u) -> Self { - Self::new(v.x().into(), v.y().into(), v.z().into()) + Self::new(v.x() as _, v.y() as _, v.z() as _) } } diff --git a/src/geom/vecu.rs b/src/geom/vecu.rs index c1da694..b55760b 100644 --- a/src/geom/vecu.rs +++ b/src/geom/vecu.rs @@ -1,7 +1,7 @@ use super::Vec3; use std::fmt::{self, Display}; -#[derive(Copy, Clone, Default, Eq, PartialEq)] +#[derive(Copy, Clone, Debug, Default, Eq, PartialEq)] pub struct Vec3u { y: u32, x: u32, diff --git a/src/mat.rs b/src/mat.rs index 464fa64..2f86a82 100644 --- a/src/mat.rs +++ b/src/mat.rs @@ -8,6 +8,9 @@ use std::cmp::Ordering; #[enum_dispatch] pub trait Material { + fn is_vacuum(&self) -> bool { + false + } /// Return \sigma, the electrical conductivity. /// For a vacuum, this is zero. For a perfect conductor, \inf. fn conductivity(&self) -> Vec3 { @@ -43,6 +46,9 @@ impl Static { } impl Material for Static { + fn is_vacuum(&self) -> bool { + self.conductivity == Default::default() + } fn conductivity(&self) -> Vec3 { self.conductivity } diff --git a/src/meas.rs b/src/meas.rs index 3b2e877..016d718 100644 --- a/src/meas.rs +++ b/src/meas.rs @@ -115,7 +115,7 @@ impl AbstractMeasurement for Energy { let f = state.feature_size(); #[allow(non_snake_case)] let dV = f*f*f; - let e = state.map_sum(|cell| { + let e = state.map_sum(|cell| { // E . D is perpetually 0 since we don't model D. // All potential energy is in the magnetic field. 0.5 * cell.h().dot(cell.b()) * dV diff --git a/src/render.rs b/src/render.rs index bbc602b..e5b12b8 100644 --- a/src/render.rs +++ b/src/render.rs @@ -70,9 +70,15 @@ impl<'a> RenderSteps<'a> { width = (width as f32 * stretch) as _; height = max_height; } - trace!("rendering at {}x{}", width, height); + trace!("rendering at {}x{} with z={}", width, height, z); 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(10.0, false, 2, |cell| { + cell.mat().conductivity().mag() + if cell.mat().is_vacuum() { + 0.0 + } else { + 5.0 + } + }); me.render_scalar_field(100.0, true, 0, |cell| cell.mat().m().mag()); if false { me.render_b_z_field();