colin/fdtd-coremem
1
0
Fork 0
Code Issues Pull Requests Projects Releases Wiki Activity

Troubleshoot toroid2.5d

Browse Source
This commit is contained in:
Colin 2020-10-09 20:41:38 -07:00
parent b9a48641c1
commit f2ba6fdd1b
8 changed files with 123 additions and 98 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

193
examples/toroid25d.rs
View File

@ -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::<Flt>() * 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::<Flt>() * 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();
}
}
}

2
src/driver.rs
View File

@ -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();
}

4
src/geom/units.rs
View File

@ -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 {

2
src/geom/vec.rs
View File

@ -199,7 +199,7 @@ impl From<(Real, Real, Real)> for Vec3 {
impl From<Vec3u> 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 _)
}
}

2
src/geom/vecu.rs
View File

@ -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,

6
src/mat.rs
View File

@ -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
}

2
src/meas.rs
View File

@ -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

10
src/render.rs
View File

@ -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();