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

Make the float depth compile-time configurable

Browse Source
This commit is contained in:
Colin
2020-09-13 22:38:00 -07:00
parent d4334565a7
commit f50dbf7d4e
9 changed files with 175 additions and 133 deletions
Download Patch File Download Diff File Expand all files Collapse all files

12
examples/em_reflection.rs
View File

@@ -5,17 +5,9 @@ fn main() {
let width = 201;
let mut driver = Driver::new(width, 101, 1e-3 /* feature size */);
driver.add_y4m_renderer("em_reflection.y4m");
driver.add_term_renderer();
// driver.add_term_renderer();
for inset in 0..20 {
for x in 0..width {
*driver.state.get_mut(x, inset).mat_mut() = mat::Static::conductor(0.1*(20.0 - inset as f64)).into();
}
for y in 0..100 {
*driver.state.get_mut(inset, y).mat_mut() = mat::Static::conductor(0.1*(20.0 - inset as f64)).into();
*driver.state.get_mut(width-1 - inset, y).mat_mut() = mat::Static::conductor(0.1*(20.0 - inset as f64)).into();
}
}
driver.add_boundary(20, 0.1);
for y in 75..100 {
for x in 0..width {
// from https://www.thoughtco.com/table-of-electrical-resistivity-conductivity-608499

1
examples/toroid2.rs
View File

@@ -16,6 +16,7 @@ fn main() {
//driver.set_steps_per_frame(40);
driver.add_y4m_renderer("toroid2.y4m");
// driver.add_term_renderer();
driver.add_measurement(meas::Label(format!("Conductivity: {}, Imax: {:.2e}", conductivity, peak_current)));
driver.add_measurement(meas::Current(width / 2 + (conductor_inner_rad + conductor_outer_rad) / 2, width / 2));
driver.add_measurement(meas::Current(width / 2 + conductor_inner_rad + 2, width / 2));
driver.add_measurement(meas::Magnetization(width / 2 + ferro_inner_rad + 2, width / 2));

12
src/driver.rs
View File

@@ -1,4 +1,4 @@
use crate::mat;
use crate::{flt, mat};
use crate::meas::{self, AbstractMeasurement};
use crate::render::{self, MultiRenderer, Renderer};
use crate::sim::SimState;
@@ -17,11 +17,11 @@ pub struct Driver {
}
impl Driver {
pub fn new(width: u32, height: u32, feature_size: f64) -> Self {
pub fn new(width: u32, height: u32, feature_size: flt::Pub) -> Self {
Driver {
state: SimState::new(width, height, feature_size),
steps_per_frame: 1,
measurements: vec![Box::new(meas::Time), Box::new(meas::Energy)],
measurements: vec![Box::new(meas::Time), Box::new(meas::Meta), Box::new(meas::Energy)],
..Default::default()
}
}
@@ -49,10 +49,10 @@ impl Driver {
self.add_renderer(render::ColorTermRenderer, "terminal");
}
pub fn add_boundary(&mut self, thickness: u32, base_conductivity: f64) {
pub fn add_boundary(&mut self, thickness: u32, base_conductivity: flt::Pub) {
for inset in 0..thickness {
let depth = thickness - inset;
let conductivity = base_conductivity * (depth*depth) as f64;
let conductivity = base_conductivity * (depth*depth) as flt::Pub;
for x in inset..self.state.width() - inset {
// left
*self.state.get_mut(x, inset).mat_mut() = mat::Static::conductor(conductivity).into();
@@ -76,7 +76,7 @@ impl Driver {
self.state.step();
self.time_spent_stepping += start_time.elapsed().unwrap();
if self.state.step_no() % (10*self.steps_per_frame) == 0 {
info!("fps: {:6.2}", (self.state.step_no() as f64) / self.time_spent_stepping.as_secs_f64());
info!("fps: {:6.2}", (self.state.step_no() as flt::Pub) / (self.time_spent_stepping.as_secs_f64() as flt::Pub));
}
}
}

72
src/geom.rs
View File

@@ -1,18 +1,18 @@
use crate::F64;
use crate::flt;
use std::ops::{Add, AddAssign, Mul, Neg, Sub};
fn real(f: f64) -> F64 {
F64::from_inner(f)
fn real(f: flt::Pub) -> flt::Priv {
flt::Priv::from_inner(f)
}
fn round(f: F64) -> F64 {
F64::from_inner(f.into_inner().round())
fn round(f: flt::Priv) -> flt::Priv {
flt::Priv::from_inner(f.into_inner().round())
}
#[derive(Copy, Clone, Debug, Default)]
pub struct Point {
pub x: F64,
pub y: F64,
pub x: flt::Priv,
pub y: flt::Priv,
}
impl Add for Point {
@@ -48,9 +48,9 @@ impl Sub for Point {
}
}
impl Mul<f64> for Point {
impl Mul<flt::Pub> for Point {
type Output = Point;
fn mul(self, s: f64) -> Point {
fn mul(self, s: flt::Pub) -> Point {
Point {
x: self.x * s,
y: self.y * s,
@@ -58,25 +58,25 @@ impl Mul<f64> for Point {
}
}
impl Into<(f64, f64)> for Point {
fn into(self) -> (f64, f64) {
impl Into<(flt::Pub, flt::Pub)> for Point {
fn into(self) -> (flt::Pub, flt::Pub) {
(self.x(), self.y())
}
}
impl Point {
pub fn new(x: f64, y: f64) -> Self {
pub fn new(x: flt::Pub, y: flt::Pub) -> Self {
Self {
x: x.into(),
y: y.into(),
}
}
pub fn x(&self) -> f64 {
pub fn x(&self) -> flt::Pub {
self.x.into()
}
pub fn y(&self) -> f64 {
pub fn y(&self) -> flt::Pub {
self.y.into()
}
@@ -87,20 +87,20 @@ impl Point {
}
}
pub fn mag_sq(&self) -> f64 {
pub fn mag_sq(&self) -> flt::Pub {
(self.x*self.x + self.y*self.y).into()
}
pub fn mag(&self) -> f64 {
pub fn mag(&self) -> flt::Pub {
self.mag_sq().sqrt()
}
pub fn with_mag(&self, new_mag: f64) -> Point {
pub fn with_mag(&self, new_mag: flt::Pub) -> Point {
if new_mag == 0.0 {
// avoid div-by-zero if self.mag() == 0 and new_mag == 0
Point::new(0.0, 0.0)
} else {
let scale = F64::from_inner(new_mag) / self.mag();
let scale = flt::Priv::from_inner(new_mag) / self.mag();
self.clone() * scale.into_inner()
}
}
@@ -130,17 +130,17 @@ impl Line {
}
}
pub fn x(&self, y: f64) -> f64 {
pub fn x(&self, y: flt::Pub) -> flt::Pub {
self.shifted(Point::new(0.0, -y)).x_intercept()
}
pub fn y(&self, x: f64) -> f64 {
pub fn y(&self, x: flt::Pub) -> flt::Pub {
(self.from.y + (real(x) - self.from.x) * self.slope()).into()
}
pub fn at_x(&self, x: f64) -> Point {
pub fn at_x(&self, x: flt::Pub) -> Point {
Point::new(x, self.y(x))
}
pub fn x_intercept(&self) -> f64 {
pub fn x_intercept(&self) -> flt::Pub {
let t = (-self.from.y) / (self.to.y - self.from.y);
(self.from.x + t * (self.to.x - self.from.x)).into()
}
@@ -149,12 +149,12 @@ impl Line {
self.to + (-self.from)
}
pub fn length_sq(&self) -> f64 {
pub fn length_sq(&self) -> flt::Pub {
let Point { x, y } = self.as_vector();
(x*x + y*y).into()
}
pub fn distance_sq(&self, pt: Point) -> f64 {
pub fn distance_sq(&self, pt: Point) -> flt::Pub {
// source: https://en.wikipedia.org/wiki/Distance_from_a_point_to_a_line#Line_defined_by_two_points
let d = self.as_vector();
let twice_area = d.y*pt.x - d.x*pt.y + self.to.x*self.from.y - self.to.y*self.from.x;
@@ -170,7 +170,7 @@ impl Line {
self.to.x > self.from.x
}
pub fn contains_x(&self, x: f64) -> bool {
pub fn contains_x(&self, x: flt::Pub) -> bool {
let x = real(x);
if self.is_ascending_x() {
x >= self.from.x && x < self.to.x
@@ -179,7 +179,7 @@ impl Line {
}
}
pub fn contains_y(&self, y: f64) -> bool {
pub fn contains_y(&self, y: flt::Pub) -> bool {
let y = real(y);
if self.is_ascending_x() {
y >= self.from.y && y < self.to.y
@@ -188,15 +188,15 @@ impl Line {
}
}
pub fn min_x(&self) -> f64 {
pub fn min_x(&self) -> flt::Pub {
self.from.x.min(self.to.x).into()
}
pub fn max_x(&self) -> f64 {
pub fn max_x(&self) -> flt::Pub {
self.from.x.max(self.to.x).into()
}
pub fn clamp_x(&self, x: f64) -> f64 {
pub fn clamp_x(&self, x: flt::Pub) -> flt::Pub {
match x {
v if v < self.min_x() => self.min_x(),
v if v > self.max_x() => self.max_x(),
@@ -204,18 +204,18 @@ impl Line {
}
}
pub fn clamp_by_x(&self, x: f64) -> Point {
pub fn clamp_by_x(&self, x: flt::Pub) -> Point {
self.at_x(self.clamp_x(x))
}
pub fn slope(&self) -> f64 {
pub fn slope(&self) -> flt::Pub {
let s = (self.to.y - self.from.y) / (self.to.x - self.from.x);
s.into()
}
/// Move the coordinate along the line toward some desired `y(x)` value.
/// The returned value could be OOB: check `self.contains_x` to find out.
pub fn move_toward_y_unclamped(&self, start_x: f64, target_y: f64) -> f64 {
pub fn move_toward_y_unclamped(&self, start_x: flt::Pub, target_y: flt::Pub) -> flt::Pub {
let start_y = self.y(start_x);
let delta_x = real(target_y - start_y) / self.slope();
(real(start_x) + delta_x).into()
@@ -249,19 +249,19 @@ impl Polygon {
})
}
pub fn min_x(&self) -> f64 {
pub fn min_x(&self) -> flt::Pub {
self.points.iter().map(|p| p.x).min().unwrap().into()
}
pub fn max_x(&self) -> f64 {
pub fn max_x(&self) -> flt::Pub {
self.points.iter().map(|p| p.x).max().unwrap().into()
}
pub fn min_y(&self) -> f64 {
pub fn min_y(&self) -> flt::Pub {
self.points.iter().map(|p| p.y).min().unwrap().into()
}
pub fn max_y(&self) -> f64 {
pub fn max_y(&self) -> flt::Pub {
self.points.iter().map(|p| p.y).max().unwrap().into()
}
}

50
src/lib.rs
View File

@@ -25,42 +25,60 @@ pub use sim::*;
// mu_0 = vacuum permeability = 1.25663706212(19)×10−6 H/m
// c_0 = speed of light in vacuum = 299792458
#[cfg(debug)]
pub(crate) type F64 = decorum::Real<f64>;
#[cfg(not(debug))]
pub(crate) type F64 = decorum::Total<f64>;
#[allow(unused)]
pub(crate) mod flt {
use decorum::{R32, R64};
pub(crate) type T64 = decorum::Total<f64>;
pub(crate) type T32 = decorum::Total<f32>;
// Use checked floats for debug and fast floats for release
#[cfg(debug)]
pub(crate) mod defaults {
use super::*;
pub(crate) type Priv = R64;
pub type Pub = f64;
}
#[cfg(not(debug))]
pub(crate) mod defaults {
use super::*;
pub(crate) type Priv = T32;
pub type Pub = f32;
}
pub use defaults::*;
}
#[allow(non_snake_case)]
pub mod consts {
use super::*;
/// Speed of light in a vacuum; m/s.
/// Also equal to 1/sqrt(epsilon_0 mu_0)
pub const C: f64 = 299792458.0;
// pub const Z0: f64 = 376.73031366857f32;
pub const C: flt::Pub = 299792458.0;
// pub const Z0: flt::Pub = 376.73031366857f32;
/// Vacuum Permeability
pub const MU0: f64 = 1.2566370621219e-6; // H/m
pub const MU0: flt::Pub = 1.2566370621219e-6; // H/m
// Vacuum Permittivity
pub const EPS0: f64 = 8.854187812813e-12; // F⋅m−1
pub const EPS0: flt::Pub = 8.854187812813e-12; // F⋅m−1
pub mod real {
use crate::F64;
pub fn C() -> F64 {
use super::*;
pub(crate) fn C() -> flt::Priv {
super::C.into()
}
pub fn C2() -> F64 {
pub(crate) fn C2() -> flt::Priv {
C() * C()
}
pub fn EPS0() -> F64 {
pub(crate) fn EPS0() -> flt::Priv {
super::EPS0.into()
}
pub fn MU0() -> F64 {
pub(crate) fn MU0() -> flt::Priv {
super::MU0.into()
}
pub fn ONE() -> F64 {
pub(crate) fn ONE() -> flt::Priv {
1.0.into()
}
pub fn TWO() -> F64 {
pub(crate) fn TWO() -> flt::Priv {
2.0.into()
}
pub fn HALF() -> F64 {
pub(crate) fn HALF() -> flt::Priv {
0.5.into()
}
}

55
src/mat.rs
View File

@@ -1,4 +1,4 @@
use crate::{CellState, consts};
use crate::{CellState, consts, flt};
use crate::geom::{Line, Point, Polygon};
use enum_dispatch::enum_dispatch;
use std::cmp::Ordering;
@@ -7,15 +7,15 @@ 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) -> f64 {
fn conductivity(&self) -> flt::Pub {
0.0
}
/// Return the magnetization in the z direction (M_z).
fn mz(&self) -> f64 {
fn mz(&self) -> flt::Pub {
0.0
}
/// Return the new h_z, and optionally change any internal state (e.g. magnetization).
fn step_h(&mut self, context: &CellState, delta_bz: f64) -> f64 {
fn step_h(&mut self, context: &CellState, delta_bz: flt::Pub) -> flt::Pub {
// B = mu0*(H+M), and in a default material M=0.
context.hz() + delta_bz / consts::MU0
}
@@ -23,12 +23,12 @@ pub trait Material {
#[derive(Clone, Default)]
pub struct Static {
pub conductivity: f64,
pub mz: f64,
pub conductivity: flt::Pub,
pub mz: flt::Pub,
}
impl Static {
pub fn conductor(conductivity: f64) -> Self {
pub fn conductor(conductivity: flt::Pub) -> Self {
Self {
conductivity,
..Default::default()
@@ -37,10 +37,10 @@ impl Static {
}
impl Material for Static {
fn conductivity(&self) -> f64 {
fn conductivity(&self) -> flt::Pub {
self.conductivity
}
fn mz(&self) -> f64 {
fn mz(&self) -> flt::Pub {
self.mz
}
}
@@ -48,19 +48,19 @@ impl Material for Static {
#[derive(Clone, Default)]
pub struct ComsolFerromagnet {
/// Instantaneous magnetization in the z direction.
pub mz: f64,
pub mz: flt::Pub,
// Ferromagnetic parameters:
/// Magnetic saturation (symmetric, so saturates to either +Ms or -Ms).
pub ms: f64,
pub ms: flt::Pub,
/// Rate at which M changes w.r.t. H.
pub xi: f64,
pub xi: flt::Pub,
}
impl Material for ComsolFerromagnet {
fn mz(&self) -> f64 {
fn mz(&self) -> flt::Pub {
self.mz
}
fn step_h(&mut self, context: &CellState, delta_bz: f64) -> f64 {
fn step_h(&mut self, context: &CellState, delta_bz: flt::Pub) -> flt::Pub {
// COMSOL ferromagnetic model: https://www.comsol.com/blogs/accessing-external-material-models-for-magnetic-simulations/
// delta(M) = xi * delta*(H) # delta(M) = M(t+1) - M(t)
// new_M = clamp(M+delta(M), -Ms, Ms)) # i.e. clamp to saturation
@@ -113,7 +113,7 @@ impl MHCurve {
/// Moves (h, m) towards some location in the MH curve where H + M = target_hm.
/// Returns `Ok((h, m))` if complete; `Err((h, m))` if there's more work to be done (call it
/// again).
fn step_toward(&self, h: f64, m: f64, target_hm: f64) -> Result<(f64, f64), (f64, f64)> {
fn step_toward(&self, h: flt::Pub, m: flt::Pub, target_hm: flt::Pub) -> Result<(flt::Pub, flt::Pub), (flt::Pub, flt::Pub)> {
let is_ascending = match target_hm.partial_cmp(&(h + m)).unwrap() {
Ordering::Greater => true,
Ordering::Less => false,
@@ -153,7 +153,7 @@ impl MHCurve {
Err(active_segment.clamp_by_x(new_h).into())
}
}
fn move_to(&self, mut h: f64, mut m: f64, target_hm: f64) -> (f64, f64) {
fn move_to(&self, mut h: flt::Pub, mut m: flt::Pub, target_hm: flt::Pub) -> (flt::Pub, flt::Pub) {
loop {
match self.step_toward(h, m, target_hm) {
Ok((x, y)) => break (x, y),
@@ -169,7 +169,8 @@ impl MHCurve {
#[derive(Clone)]
pub struct PiecewiseLinearFerromagnet {
/// Instantaneous magnetization in the z direction.
pub mz: f64,
pub mz: flt::Pub,
pub conductivity: flt::Pub,
mh_curve: MHCurve,
}
@@ -177,23 +178,27 @@ impl PiecewiseLinearFerromagnet {
/// Construct a B(H) curve from a sweep from H = 0 to Hs and back down to H = 0.
/// The curve below B = 0 is derived by symmetry.
/// Points are (H, B).
pub fn from_bh(points: &[(f64, f64)]) -> Self {
pub fn from_bh(points: &[(flt::Pub, flt::Pub)]) -> Self {
let mh_points: Vec<Point> = points.iter().cloned().map(|(h, b)| {
Point::new(h, b / consts::MU0 - h)
}).collect();
Self {
mz: 0.0,
conductivity: 0.0,
mh_curve: MHCurve::new(&*mh_points),
}
}
}
impl Material for PiecewiseLinearFerromagnet {
fn mz(&self) -> f64 {
fn mz(&self) -> flt::Pub {
self.mz
}
fn step_h(&mut self, context: &CellState, delta_bz: f64) -> f64 {
fn conductivity(&self) -> flt::Pub {
self.conductivity
}
fn step_h(&mut self, context: &CellState, delta_bz: flt::Pub) -> flt::Pub {
let (h, m) = (context.hz(), self.mz);
let target_hm = h + m + delta_bz/consts::MU0;
@@ -223,7 +228,7 @@ pub mod db {
use super::*;
/// https://www.ferroxcube.com/upload/media/product/file/MDS/3r1.pdf
pub fn ferroxcube_3r1() -> GenericMaterial {
PiecewiseLinearFerromagnet::from_bh(&[
let mut m = PiecewiseLinearFerromagnet::from_bh(&[
( 35.0, 0.0),
( 50.0, 0.250),
( 100.0, 0.325),
@@ -234,7 +239,9 @@ pub mod db {
( 100.0, 0.345),
( 50.0, 0.340),
( 0.0, 0.325),
]).into()
]);
m.conductivity = 1e3;
m.into()
}
}
@@ -258,7 +265,7 @@ mod test {
])
}
fn assert_step_toward_symmetric(h: f64, m: f64, target_mh: f64, target: Result<(f64, f64), (f64, f64)>) {
fn assert_step_toward_symmetric(h: flt::Pub, m: flt::Pub, target_mh: flt::Pub, target: Result<(flt::Pub, flt::Pub), (flt::Pub, flt::Pub)>) {
let curve = mh_curve_for_test();
let neg_target = match target {
Ok((a, b)) => Ok((-a, -b)),
@@ -268,7 +275,7 @@ mod test {
assert_eq!(curve.step_toward(-h, -m, -target_mh), neg_target);
}
fn assert_move_to_symmetric(h: f64, m: f64, target_mh: f64, target: (f64, f64)) {
fn assert_move_to_symmetric(h: flt::Pub, m: flt::Pub, target_mh: flt::Pub, target: (flt::Pub, flt::Pub)) {
let curve = mh_curve_for_test();
assert_eq!(curve.move_to(h, m, target_mh), target);
assert_eq!(curve.move_to(-h, -m, -target_mh), (-target.0, -target.1));

24
src/meas.rs
View File

@@ -9,7 +9,29 @@ pub struct Time;
impl AbstractMeasurement for Time {
fn eval(&self, state: &SimSnapshot) -> String {
format!("time: {:.3e} (step {})", state.time(), state.step_no())
format!("{:.3e}s (step {})", state.time(), state.step_no())
}
}
pub struct Meta;
impl AbstractMeasurement for Meta {
fn eval(&self, state: &SimSnapshot) -> String {
format!("{}x{} feat: {:.1e}m", state.width(), state.height(), state.feature_size())
}
}
pub struct Label(pub String);
impl Label {
pub fn new<S: Into<String>>(s: S) -> Self {
Self(s.into())
}
}
impl AbstractMeasurement for Label {
fn eval(&self, _state: &SimSnapshot) -> String {
self.0.clone()
}
}

30
src/render.rs
View File

@@ -1,10 +1,10 @@
use ansi_term::Color::RGB;
use crate::geom::Point;
use crate::{F64, Material as _, SimSnapshot, SimState};
use crate::{flt, Material as _, SimSnapshot, SimState};
use crate::mat;
use crate::sim::Cell;
use crate::meas::AbstractMeasurement;
use font8x8::{BASIC_FONTS, UnicodeFonts as _};
use font8x8::{BASIC_FONTS, GREEK_FONTS, UnicodeFonts as _};
use image::{RgbImage, Rgb};
use imageproc::{pixelops, drawing};
use std::fs::File;
@@ -14,7 +14,7 @@ use y4m;
/// Accept a value from (-\inf, \inf) and return a value in (-1, 1).
/// If the input is equal to `typical`, it will be mapped to 0.5.
/// If the input is equal to -`typical`, it will be mapped to -0.5.
fn scale_signed(x: f64, typical: f64) -> f64 {
fn scale_signed(x: flt::Pub, typical: flt::Pub) -> flt::Pub {
if x >= 0.0 {
scale_unsigned(x, typical)
} else {
@@ -24,29 +24,29 @@ fn scale_signed(x: f64, typical: f64) -> f64 {
/// Accept a value from [0, \inf) and return a value in [0, 1).
/// If the input is equal to `typical`, it will be mapped to 0.5.
fn scale_unsigned(x: f64, typical: f64) -> f64 {
fn scale_unsigned(x: flt::Pub, typical: flt::Pub) -> flt::Pub {
// f(0) => 0
// f(1) => 0.5
// f(\inf) => 1
// f(x) = 1 - 1/(x+1)
let x = F64::from_inner(x);
let typ = F64::from_inner(typical);
let y = F64::from_inner(1.0) - F64::from_inner(1.0)/(x/typ + 1.0);
let x = flt::Priv::from_inner(x);
let typ = flt::Priv::from_inner(typical);
let y = flt::Priv::from_inner(1.0) - flt::Priv::from_inner(1.0)/(x/typ + 1.0);
y.into()
}
fn scale_signed_to_u8(x: f64, typ: f64) -> u8 {
fn scale_signed_to_u8(x: flt::Pub, typ: flt::Pub) -> u8 {
let norm = 128.0 + 128.0*scale_signed(x, typ);
norm as _
}
fn scale_unsigned_to_u8(x: f64, typ: f64) -> u8 {
fn scale_unsigned_to_u8(x: flt::Pub, typ: flt::Pub) -> u8 {
let norm = 256.0*scale_unsigned(x, typ);
norm as _
}
/// Scale a vector to have magnitude between [0, 1).
fn scale_vector(x: Point, typical_mag: f64) -> Point {
fn scale_vector(x: Point, typical_mag: flt::Pub) -> Point {
let new_mag = scale_unsigned(x.mag(), typical_mag);
x.with_mag(new_mag)
}
@@ -88,7 +88,7 @@ impl<'a> RenderSteps<'a> {
}
}
}
fn render_vector_field<F: Fn(&Cell<mat::Static>) -> Point>(&mut self, color: Rgb<u8>, typical: f64, measure: F) {
fn render_vector_field<F: Fn(&Cell<mat::Static>) -> Point>(&mut self, color: Rgb<u8>, typical: flt::Pub, measure: F) {
let w = self.sim.width();
let h = self.sim.height();
let vec_spacing = 10;
@@ -98,7 +98,7 @@ impl<'a> RenderSteps<'a> {
let vec = self.field_vector(x, y, vec_spacing, &measure);
let norm_vec = scale_vector(vec, typical);
let alpha = 0.7*scale_unsigned(vec.mag_sq(), typical * 5.0);
let vec = norm_vec * (vec_spacing as f64);
let vec = norm_vec * (vec_spacing as flt::Pub);
let center = Point::new(x as _, y as _) + Point::new(vec_spacing as _, vec_spacing as _)*0.5;
self.im.draw_field_arrow(center, vec, color, alpha as f32);
}
@@ -120,7 +120,9 @@ impl<'a> RenderSteps<'a> {
for (meas_no, m) in self.meas.iter().enumerate() {
let meas_string = m.eval(self.sim);
for (i, c) in meas_string.chars().enumerate() {
let glyph = BASIC_FONTS.get(c).unwrap();
let glyph = BASIC_FONTS.get(c)
.or_else(|| GREEK_FONTS.get(c))
.unwrap_or_else(|| BASIC_FONTS.get('?').unwrap());
for (y, bmp) in glyph.iter().enumerate() {
for x in 0..8 {
if (bmp & 1 << x) != 0 {
@@ -155,7 +157,7 @@ impl<'a> RenderSteps<'a> {
// avoid division by zero
Point::new(0.0, 0.0)
} else {
field * (1.0 / ((xw*yw) as f64))
field * (1.0 / ((xw*yw) as flt::Pub))
}
}
}

52
src/sim.rs
View File

@@ -1,4 +1,4 @@
use crate::{F64, consts};
use crate::{flt, consts};
use crate::geom::Point;
use crate::mat::{self, GenericMaterial, Material};
@@ -12,12 +12,12 @@ pub type SimSnapshot = SimState<mat::Static>;
pub struct SimState<M=GenericMaterial> {
cells: Array2<Cell<M>>,
scratch: Array2<Cell<M>>,
feature_size: F64,
feature_size: flt::Priv,
step_no: u64,
}
impl<M: Material + Default> SimState<M> {
pub fn new(width: u32, height: u32, feature_size: f64) -> Self {
pub fn new(width: u32, height: u32, feature_size: flt::Pub) -> Self {
Self {
cells: Array2::default((height as _, width as _)),
scratch: Array2::default((height as _, width as _)),
@@ -64,7 +64,7 @@ impl<M: Material + Default + Send + Sync> SimState<M> {
pub fn snapshot(&self, dec: u32) -> SimSnapshot {
let rows = self.height() / dec;
let cols = self.width() / dec;
let sample_area = (dec * dec) as f64;
let sample_area = (dec * dec) as flt::Pub;
// XXX Zip::par_apply_collect requires the Cell type to be Copy,
// so we assign into a pre-allocated type instead.
@@ -102,28 +102,28 @@ impl<M: Material + Default + Send + Sync> SimState<M> {
}
impl<M: Material> SimState<M> {
pub fn impulse_bz(&mut self, x: u32, y: u32, bz: f64) {
pub fn impulse_bz(&mut self, x: u32, y: u32, bz: flt::Pub) {
self.get_mut(x, y).impulse_bz(bz);
}
}
impl<M> SimState<M> {
pub fn time(&self) -> f64 {
(self.timestep() * self.step_no as f64).into()
pub fn time(&self) -> flt::Pub {
(self.timestep() * self.step_no as flt::Pub).into()
}
pub fn step_no(&self) -> u64 {
self.step_no
}
pub fn feature_size(&self) -> f64 {
pub fn feature_size(&self) -> flt::Pub {
self.feature_size.into()
}
pub fn impulse_ex(&mut self, x: u32, y: u32, ex: f64) {
pub fn impulse_ex(&mut self, x: u32, y: u32, ex: flt::Pub) {
self.get_mut(x, y).state.ex += ex;
}
pub fn impulse_ey(&mut self, x: u32, y: u32, ey: f64) {
pub fn impulse_ey(&mut self, x: u32, y: u32, ey: flt::Pub) {
self.get_mut(x, y).state.ey += ey;
}
@@ -140,7 +140,7 @@ impl<M> SimState<M> {
&mut self.cells[[y as _, x as _]]
}
fn timestep(&self) -> F64 {
fn timestep(&self) -> flt::Priv {
self.feature_size / consts::real::C()
}
}
@@ -170,16 +170,16 @@ pub struct Cell<M> {
}
impl<M> Cell<M> {
pub fn ex(&self) -> f64 {
pub fn ex(&self) -> flt::Pub {
self.state.ex()
}
pub fn ey(&self) -> f64 {
pub fn ey(&self) -> flt::Pub {
self.state.ey()
}
pub fn e(&self) -> Point {
Point::new(self.ex(), self.ey())
}
pub fn hz(&self) -> f64 {
pub fn hz(&self) -> flt::Pub {
self.state.hz()
}
pub fn mat(&self) -> &M {
@@ -191,7 +191,7 @@ impl<M> Cell<M> {
}
impl<M: Material> Cell<M> {
pub fn bz(&self) -> f64 {
pub fn bz(&self) -> flt::Pub {
consts::MU0 * (self.hz() + self.mat.mz())
}
@@ -200,16 +200,16 @@ impl<M: Material> Cell<M> {
Point::new(self.ex()*conductivity, self.ey()*conductivity)
}
fn bz_to_hz(&self, bz: f64) -> f64 {
fn bz_to_hz(&self, bz: flt::Pub) -> flt::Pub {
// B = mu0*(H + M) => H = B/mu0 - M
bz/consts::MU0 - self.mat.mz()
}
fn impulse_bz(&mut self, delta_bz: f64) {
fn impulse_bz(&mut self, delta_bz: flt::Pub) {
self.state.hz = self.mat.step_h(&self.state, delta_bz).into();
}
fn step_h(mut self, right: &Self, down: &Self, _delta_t: F64) -> Self {
fn step_h(mut self, right: &Self, down: &Self, _delta_t: flt::Priv) -> Self {
// ```tex
// Maxwell-Faraday equation: $\nabla x E = -dB/dt$
// Expand: $dE_y/dx - dE_x/dy = -dB_z/dt$
@@ -235,7 +235,7 @@ impl<M: Material> Cell<M> {
/// delta_t = timestep covered by this function. i.e. it should be half the timestep of the simulation
/// feature_size = how many units apart is the center of each adjacent cell on the grid.
fn step_e(self, left: &Self, up: &Self, delta_t: F64, feature_size: F64) -> Self {
fn step_e(self, left: &Self, up: &Self, delta_t: flt::Priv, feature_size: flt::Priv) -> Self {
// ```tex
// Ampere's circuital law with Maxwell's addition, in SI units ("macroscopic version"):
// $\nabla x H = J_f + dD/dt$ where $J_f$ = current density = $\sigma E$, $\sigma$ being a material parameter ("conductivity")
@@ -262,7 +262,7 @@ impl<M: Material> Cell<M> {
use consts::real::{EPS0, ONE, TWO};
let sigma: F64 = self.mat.conductivity().into();
let sigma: flt::Priv = self.mat.conductivity().into();
// XXX not obvious that bz_to_hz is sensible
// let delta_hz_y = self.hz() - self.bz_to_hz(up.bz());
@@ -288,19 +288,19 @@ impl<M: Material> Cell<M> {
#[derive(Copy, Clone, Default)]
pub struct CellState {
ex: F64,
ey: F64,
hz: F64,
ex: flt::Priv,
ey: flt::Priv,
hz: flt::Priv,
}
impl CellState {
pub fn ex(&self) -> f64 {
pub fn ex(&self) -> flt::Pub {
self.ex.into()
}
pub fn ey(&self) -> f64 {
pub fn ey(&self) -> flt::Pub {
self.ey.into()
}
pub fn hz(&self) -> f64 {
pub fn hz(&self) -> flt::Pub {
self.hz.into()
}
}