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

Add a 'Measurement' concept, e.g. to display current at a point.

Browse Source
This commit is contained in:
Colin
2020-09-07 20:14:41 -07:00
parent fbca8622df
commit d611e822dd
7 changed files with 138 additions and 47 deletions
Download Patch File Download Diff File Expand all files Collapse all files

1
Cargo.toml
View File

@@ -10,6 +10,7 @@ edition = "2018"
ansi_term = "0.12"
decorum = "0.3"
enum_dispatch = "0.3"
font8x8 = "0.2"
image = "0.23"
imageproc = "0.21"
ndarray = "0.13"

58
examples/ferromagnet.rs
View File

@@ -1,16 +1,21 @@
use coremem::{Driver, mat};
use coremem::{Driver, mat, meas};
fn main() {
let width = 1300;
let height = 1000;
let mut driver = Driver::new(width, height, 1e-5 /* feature size */);
driver.set_steps_per_frame(40);
let width = 500;
let height = 500;
let feat_size = 1e-5; // feature size
let conductivity = 1.0e3;
let mut driver = Driver::new(width, height, feat_size);
driver.set_steps_per_frame(8);
//driver.set_steps_per_frame(40);
driver.add_y4m_renderer("ferromagnet.y4m");
driver.add_term_renderer();
// driver.add_term_renderer();
driver.add_measurement(meas::Current(225, 250));
driver.add_measurement(meas::Current(300, 250));
for y in 0..height {
for x in 500..600 {
*driver.state.get_mut(x, y).mat_mut() = mat::Static::conductor(1.0e1).into();
for x in 200..250 {
*driver.state.get_mut(x, y).mat_mut() = mat::Static::conductor(conductivity).into();
}
// for x in 30..40 {
// *state.get_mut(x, y).mat_mut() = mat::Conductor { conductivity: 1.0e8 }.into();
@@ -20,12 +25,12 @@ fn main() {
// *state.get_mut(x, y).mat_mut() = mat::Conductor { conductivity: 1.0e8 }.into();
// }
// }
for x in 720..800 {
*driver.state.get_mut(x, y).mat_mut() = mat::Static::conductor(1.0e1).into();
for x in 280..330 {
*driver.state.get_mut(x, y).mat_mut() = mat::Static::conductor(conductivity).into();
}
}
for y in 400..600 {
for x in 620..700 {
for y in 200..300 {
for x in 260..270 {
*driver.state.get_mut(x, y).mat_mut() = mat::PiecewiseLinearFerromagnet::from_bh(&[
( 35.0, 0.0),
( 50.0, 0.250),
@@ -53,18 +58,33 @@ fn main() {
// 400..=440 => -1e6,
// _ => 0.0
//};
let imp = if driver.state.step_no() < 50 {
250000.0 * ((driver.state.step_no() as f64)*0.02*std::f64::consts::PI).sin()
} else {
0.0
// let v = if driver.state.step_no() < 50 {
// 2.5 * ((driver.state.step_no() as f64)*0.02*std::f64::consts::PI).sin()
// } else {
// 0.0
// };
let drive_current = match driver.state.step_no() {
0..=1000 => 2.5e7,
2000..=3000 => -2.5e7,
_ => 0.0,
};
// state.impulse_ex(50, 50, imp);
// 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 100..height-100 {
for x in 520..580 {
driver.state.impulse_ey(x, y, imp as _);
// for y in 100..height-100 {
// for x in 200..250 {
// driver.state.impulse_ey(x, y, imp as _);
// }
// }
// E = V/M
//let e = v/(2.0*feat_size);
let e = drive_current/conductivity;
for x in 200..250 {
for y in 100..height-100 {
// driver.state.impulse_ey(x, 102, e);
// driver.state.impulse_ey(x, height - 103, e);
driver.state.impulse_ey(x, y, e);
}
}
driver.step();

12
src/driver.rs
View File

@@ -1,4 +1,5 @@
use crate::mat;
use crate::meas::{self, AbstractMeasurement};
use crate::render::{self, MultiRenderer, Renderer};
use crate::sim::SimState;
@@ -11,6 +12,7 @@ pub struct Driver {
renderer: MultiRenderer,
steps_per_frame: u64,
time_spent_stepping: Duration,
measurements: Vec<Box<dyn AbstractMeasurement>>,
}
impl Driver {
@@ -18,10 +20,15 @@ impl Driver {
Driver {
state: SimState::new(width, height, feature_size),
steps_per_frame: 1,
measurements: vec![Box::new(meas::Time)],
..Default::default()
}
}
pub fn add_measurement<M: AbstractMeasurement + 'static>(&mut self, m: M) {
self.measurements.push(Box::new(m));
}
pub fn set_steps_per_frame(&mut self, steps_per_frame: u64) {
self.steps_per_frame = steps_per_frame;
}
@@ -42,7 +49,7 @@ impl Driver {
for inset in 0..thickness {
let depth = thickness - inset;
// TODO: tune a scalar multiplier on this value
let conductivity = (depth*depth) as f64;
let conductivity = 0.1 * (depth*depth) as f64;
for x in inset..self.state.width() - inset {
// left
*self.state.get_mut(x, inset).mat_mut() = mat::Static::conductor(conductivity).into();
@@ -61,10 +68,11 @@ impl Driver {
pub fn step(&mut self) {
let start_time = SystemTime::now();
if self.state.step_no() % self.steps_per_frame == 0 {
self.renderer.render(&self.state);
self.renderer.render(&self.state, &*self.measurements);
}
self.state.step();
self.time_spent_stepping += start_time.elapsed().unwrap();
println!("fps: {}", (self.state.step_no() as f64) / self.time_spent_stepping.as_secs_f64());
}
}

1
src/lib.rs
View File

@@ -10,6 +10,7 @@ use decorum::R64;
pub mod driver;
pub mod geom;
pub mod mat;
pub mod meas;
pub mod render;
pub mod sim;

22
src/meas.rs Normal file
View File

@@ -0,0 +1,22 @@
use crate::sim::SimSnapshot;
pub trait AbstractMeasurement {
fn eval(&self, state: &SimSnapshot) -> String;
}
pub struct Time;
impl AbstractMeasurement for Time {
fn eval(&self, state: &SimSnapshot) -> String {
format!("time: {:.3e} (step {})", state.time(), state.step_no())
}
}
pub struct Current(pub u32, pub u32);
impl AbstractMeasurement for Current {
fn eval(&self, state: &SimSnapshot) -> String {
let current = state.get(self.0, self.1).current();
format!("current({}, {}): ({:.2e}, {:.2e})", self.0, self.1, current.x, current.y)
}
}

86
src/render.rs
View File

@@ -1,7 +1,9 @@
use ansi_term::Color::RGB;
use crate::geom::Point;
use crate::{Material as _, SimSnapshot, SimState};
use crate::meas::AbstractMeasurement;
use decorum::R64;
use font8x8::{BASIC_FONTS, UnicodeFonts as _};
use image::{RgbImage, Rgb};
use imageproc::{pixelops, drawing};
use std::fs::File;
@@ -48,43 +50,75 @@ fn scale_vector(x: Point, typical_mag: f64) -> Point {
x.with_mag(new_mag)
}
trait SimSnapshotRenderExt {
fn to_image(&self) -> RgbImage;
fn e_vector(&self, xidx: u32, yidx: u32, size: u32) -> Point;
}
struct RenderSteps(RgbImage);
impl SimSnapshotRenderExt for SimSnapshot {
fn to_image(&self) -> RgbImage {
impl RenderSteps {
fn render(state: &SimSnapshot, measurements: &[Box<dyn AbstractMeasurement>]) -> RgbImage {
let w = state.width();
let h = state.height();
let mut me = RenderSteps(RgbImage::new(w, h));
me.render_b_field(state);
me.render_e_field(state);
me.render_measurements(state, measurements);
me.0
}
fn width(&self) -> u32 {
self.0.width()
}
fn height(&self) -> u32 {
self.0.height()
}
fn render_b_field(&mut self, state: &SimSnapshot) {
let w = self.width();
let h = self.height();
let evec_spacing = 10;
let mut image = RgbImage::new(w, h);
for y in 0..h {
for x in 0..w {
let cell = self.get(x, y);
let cell = state.get(x, y);
let r = scale_signed_to_u8(cell.mat().mz(), 100.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]));
self.0.put_pixel(x, y, Rgb([r, g, b]));
}
}
}
fn render_e_field(&mut self, state: &SimSnapshot) {
let w = self.width();
let h = self.height();
let evec_spacing = 10;
for y in 0..h {
for x in 0..w {
if x % evec_spacing == 0 && y % evec_spacing == 0 {
let evec = self.e_vector(x, y, evec_spacing);
let evec = self.e_vector(state, x, y, evec_spacing);
let norm_vec = scale_vector(evec, 100.0);
let alpha = scale_unsigned(evec.mag_sq(), 500.0);
let vec = norm_vec * (evec_spacing as f64);
let center = Point::new(x as _, y as _) + Point::new(evec_spacing as _, evec_spacing as _)*0.5;
image.draw_field_arrow(center, vec, Rgb([0xff, 0xff, 0xff]), alpha as f32);
self.0.draw_field_arrow(center, vec, Rgb([0xff, 0xff, 0xff]), alpha as f32);
}
}
}
}
fn render_measurements(&mut self, state: &SimSnapshot, measurements: &[Box<dyn AbstractMeasurement>]) {
for (meas_no, m) in measurements.iter().enumerate() {
let meas_string = m.eval(state);
for (i, c) in meas_string.chars().enumerate() {
let glyph = BASIC_FONTS.get(c).unwrap();
for (y, bmp) in glyph.iter().enumerate() {
for x in 0..8 {
if (bmp & 1 << x) != 0 {
let real_x = 2 + i as u32*8 + x;
let real_y = y as u32 + self.height() - 10 - meas_no as u32 * 8;
if real_x < self.width() {
self.0.put_pixel(real_x, real_y, Rgb([0, 0, 0]));
}
}
}
}
}
}
image
}
fn e_vector(&self, xidx: u32, yidx: u32, size: u32) -> Point {
fn e_vector(&self, state: &SimSnapshot, xidx: u32, yidx: u32, size: u32) -> Point {
let mut e = Point::default();
let w = self.width();
let h = self.height();
@@ -94,7 +128,7 @@ impl SimSnapshotRenderExt for SimSnapshot {
let yend = (ystart + size).min(h);
for y in ystart..yend {
for x in xstart..xend {
e += self.get(x, y).e();
e += state.get(x, y).e();
}
}
let xw = xend - xstart;
@@ -127,18 +161,18 @@ impl ImageRenderExt for RgbImage {
}
pub trait Renderer {
fn render(&mut self, state: &SimSnapshot) {
self.render_with_image(state, &state.to_image());
fn render(&mut self, state: &SimSnapshot, measurements: &[Box<dyn AbstractMeasurement>]) {
self.render_with_image(state, &RenderSteps::render(state, measurements));
}
fn render_with_image(&mut self, state: &SimSnapshot, _im: &RgbImage) {
self.render(state);
self.render(state, &[]);
}
}
pub struct NumericTermRenderer;
impl Renderer for NumericTermRenderer {
fn render(&mut self, state: &SimSnapshot) {
fn render(&mut self, state: &SimSnapshot, _measurements: &[Box<dyn AbstractMeasurement>]) {
for y in 0..state.height() {
for x in 0..state.width() {
let cell = state.get(x, y);
@@ -236,19 +270,19 @@ impl MultiRenderer {
self
}
pub fn render(&mut self, state: &SimState) {
//let max_width = 1980; //< TODO: make configurable
let max_width = 200;
pub fn render(&mut self, state: &SimState, measurements: &[Box<dyn AbstractMeasurement>]) {
let max_width = 1980; //< TODO: make configurable
// let max_width = 200;
let dec = (state.width() + max_width - 1) / max_width;
let snap = state.snapshot(dec);
Renderer::render(self, &snap);
Renderer::render(self, &snap, measurements);
}
}
impl Renderer for MultiRenderer {
fn render(&mut self, state: &SimSnapshot) {
fn render(&mut self, state: &SimSnapshot, measurements: &[Box<dyn AbstractMeasurement>]) {
if self.renderers.len() != 0 {
self.render_with_image(state, &state.to_image());
self.render_with_image(state, &RenderSteps::render(state, measurements));
}
}

5
src/sim.rs
View File

@@ -211,6 +211,11 @@ impl<M: Material> Cell<M> {
consts::MU0 * (self.hz() + self.mat.mz())
}
pub fn current(&self) -> Point {
let conductivity = self.mat.conductivity();
Point::new(self.ex()*conductivity, self.ey()*conductivity)
}
fn bz_to_hz(&self, bz: f64) -> f64 {
// B = mu0*(H + M) => H = B/mu0 - M
bz/consts::MU0 - self.mat.mz()