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

Add units to the different Vec3 "types" (Meters and Index, so far)

Browse Source
This commit is contained in:
Colin
2020-09-26 17:05:15 -07:00
parent 00f22c95d0
commit 55d6171325
9 changed files with 180 additions and 107 deletions
Download Patch File Download Diff File Expand all files Collapse all files

30
examples/toroid25d.rs
View File

@@ -1,5 +1,5 @@
use coremem::{consts, Driver, Flt, mat, meas};
use coremem::geom::{Coord, CylinderZ, Vec2, Vec3};
use coremem::{Driver, Flt, mat, meas};
use coremem::geom::{Coord, CylinderZ, Index, Meters, Vec2, Vec3};
use coremem::stim::{Stimulus, Sinusoid};
fn main() {
@@ -22,14 +22,14 @@ fn main() {
let width_px = from_m(width);
let depth_px = from_m(depth);
let size_px = (width_px, width_px, depth_px).into();
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-{}A--radii{}um-{}um-{}um.y4m",
std::mem::size_of::<Flt>() * 8,
size_px,
*size_px,
m_to_um(feat_size),
peak_current as u32,
m_to_um(conductor_outer_rad),
@@ -43,31 +43,31 @@ fn main() {
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(
(half_width + ferro_inner_rad + 2.0*feat_size, half_width, half_depth).into()
Meters((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, half_depth).into()
Meters((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, half_depth).into()
Meters((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, half_depth).into()
Meters((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, half_depth).into()
Meters((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, half_depth).into()
Meters((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, half_depth).into()
Meters((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, half_depth).into()
Meters((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, half_depth).into()
Meters((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 _);
@@ -75,7 +75,7 @@ fn main() {
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, z_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) {
@@ -90,7 +90,7 @@ fn main() {
}
let boundary_xy = half_width - ferro_outer_rad - buffer;
println!("boundary: {}um", m_to_um(boundary_xy));
let boundary = Coord::new(from_m(boundary_xy), from_m(boundary_xy), 20);
let boundary = Index((from_m(boundary_xy), from_m(boundary_xy), 20).into());
driver.add_upml_boundary(boundary);
driver.add_stimulus(Stimulus::new(

11
src/driver.rs
View File

@@ -1,6 +1,6 @@
use crate::{flt::Flt, mat};
use crate::consts;
use crate::geom::{Coord, Vec3};
use crate::geom::{Coord, Index, Vec3, Vec3u};
use crate::meas::{self, AbstractMeasurement};
use crate::render::{self, MultiRenderer, Renderer};
use crate::sim::{GenericSim as _, SimState};
@@ -23,9 +23,9 @@ pub struct Driver {
}
impl Driver {
pub fn new(size: Coord, feature_size: Flt) -> Self {
pub fn new<C: Coord>(size: C, feature_size: Flt) -> Self {
Driver {
state: SimState::new(size, feature_size),
state: SimState::new(size.to_index(feature_size), feature_size),
renderer: Default::default(),
steps_per_frame: 1,
time_spent_stepping: Default::default(),
@@ -83,8 +83,9 @@ impl Driver {
// }
// }
pub fn add_upml_boundary(&mut self, thickness: Coord) {
pub fn add_upml_boundary<C: Coord>(&mut self, thickness: C) {
// Based on explanation here (slide 63): https://empossible.net/wp-content/uploads/2020/01/Lecture-The-Perfectly-Matched-Layer.pdf
let thickness = thickness.to_index(self.state.feature_size());
let h = self.state.height();
let w = self.state.width();
let d = self.state.depth();
@@ -119,7 +120,7 @@ impl Driver {
let cond_y = scale * (depth_y as Flt/thickness.y() as Flt).powf(3.0);
let cond_z = scale * (depth_z as Flt/thickness.z() as Flt).powf(3.0);
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();
*self.state.get_mut(Index(Vec3u::new(x, y, z))).mat_mut() = conductor.into();
}
}
}

11
src/geom/mod.rs
View File

@@ -1,8 +1,10 @@
mod coord;
mod units;
mod vec;
mod vecu;
pub use coord::Coord;
pub use units::{Coord, Meters, Index};
pub use vec::{Vec2, Vec3};
pub use vecu::Vec3u;
use crate::flt::{Flt, Real};
use std::fmt::{self, Display};
@@ -180,7 +182,7 @@ impl Polygon2d {
}
pub trait Region {
fn contains(&self, p: Vec3) -> bool;
fn contains(&self, p: Meters) -> bool;
}
#[derive(Copy, Clone)]
@@ -190,6 +192,7 @@ pub struct CylinderZ {
}
impl CylinderZ {
// TODO: should be unit-typed
pub fn new(center: Vec2, radius: Flt) -> Self {
Self {
center,
@@ -199,7 +202,7 @@ impl CylinderZ {
}
impl Region for CylinderZ {
fn contains(&self, p: Vec3) -> bool {
fn contains(&self, p: Meters) -> bool {
p.xy().distance_sq(self.center) <= (self.radius * self.radius).into()
}
}

60
src/geom/units.rs Normal file
View File

@@ -0,0 +1,60 @@
use crate::flt::Flt;
use super::{Vec3, Vec3u};
use std::ops::Deref;
pub trait Coord: Copy + Clone {
fn to_meters(&self, feature_size: Flt) -> Meters;
fn to_index(&self, feature_size: Flt) -> Index;
fn from_meters(other: Meters, feature_size: Flt) -> Self;
fn from_index(other: Index, feature_size: Flt) -> Self;
}
#[derive(Copy, Clone)]
pub struct Meters(pub Vec3);
impl Coord for Meters {
fn to_meters(&self, _feature_size: Flt) -> Meters {
*self
}
fn to_index(&self, feature_size: Flt) -> Index {
Index((self.0 / feature_size).into())
}
fn from_meters(other: Meters, _feature_size: Flt) -> Self {
other
}
fn from_index(other: Index, feature_size: Flt) -> Self {
other.to_meters(feature_size)
}
}
impl Deref for Meters {
type Target = Vec3;
fn deref(&self) -> &Vec3 {
&self.0
}
}
#[derive(Copy, Clone)]
pub struct Index(pub Vec3u);
impl Coord for Index {
fn to_meters(&self, feature_size: Flt) -> Meters {
Meters(Vec3::from(self.0) * feature_size)
}
fn to_index(&self, _feature_size: Flt) -> Index {
*self
}
fn from_meters(other: Meters, feature_size: Flt) -> Self {
other.to_index(feature_size)
}
fn from_index(other: Index, _feature_size: Flt) -> Self {
other
}
}
impl Deref for Index {
type Target = Vec3u;
fn deref(&self) -> &Vec3u {
&self.0
}
}

8
src/geom/vec.rs
View File

@@ -1,5 +1,5 @@
use crate::flt::{Flt, Real};
use super::Coord;
use super::Vec3u;
use std::convert::From;
use std::iter::Sum;
use std::ops::{Add, AddAssign, Div, DivAssign, Mul, MulAssign, Neg, Sub};
@@ -197,9 +197,9 @@ impl From<(Real, Real, Real)> for Vec3 {
}
}
impl From<Coord> for Vec3 {
fn from(c: Coord) -> Self {
Self::new(c.x().into(), c.y().into(), c.z().into())
impl From<Vec3u> for Vec3 {
fn from(v: Vec3u) -> Self {
Self::new(v.x().into(), v.y().into(), v.z().into())
}
}

15
src/geom/coord.rs → src/geom/vecu.rs
View File

@@ -1,12 +1,13 @@
use super::Vec3;
use std::fmt::{self, Display};
#[derive(Copy, Clone, Default, Eq, PartialEq)]
pub struct Coord {
pub struct Vec3u {
y: u32,
x: u32,
z: u32,
}
impl Coord {
impl Vec3u {
pub fn new(x: u32, y: u32, z: u32) -> Self {
Self {
x,
@@ -29,13 +30,19 @@ impl Coord {
}
}
impl From<(u32, u32, u32)> for Coord {
impl From<(u32, u32, u32)> for Vec3u {
fn from((x, y, z): (u32, u32, u32)) -> Self {
Self::new(x, y, z)
}
}
impl Display for Coord {
impl From<Vec3> for Vec3u {
fn from(v: Vec3) -> Self {
Self::new(v.x() as _, v.y() as _, v.z() as _)
}
}
impl Display for Vec3u {
fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
write!(f, "({}, {}, {})", self.x, self.y, self.z)
}

20
src/meas.rs
View File

@@ -1,5 +1,5 @@
use crate::flt::Flt;
use crate::geom::{Coord, Vec3, Region};
use crate::geom::{Meters, Region};
use crate::mat::Material as _;
use crate::sim::{Cell, GenericSim};
use std::fmt::Display;
@@ -39,12 +39,14 @@ impl AbstractMeasurement for Label {
}
}
fn sum_over_region<T: Default + Sum<T>, R: Region, F: Fn(Coord, &Cell) -> T>(state: &dyn GenericSim, r: &R, f: F) -> T {
fn sum_over_region<T: Default + Sum<T>, R: Region, F: Fn(Meters, &Cell) -> T>(state: &dyn GenericSim, r: &R, f: F) -> T {
// TODO: z coordinate?
state.map_sum_enumerated(|coord, cell| if r.contains(Vec3::from(coord)*state.feature_size()) {
state.map_sum_enumerated(|coord, cell| {
if r.contains(coord) {
f(coord, cell)
} else {
Default::default()
}
})
}
@@ -57,13 +59,13 @@ impl<R: Region + Display + Sync> AbstractMeasurement for Current<R> {
}
}
fn loc(v: Vec3) -> String {
let (x, y, z): (Flt, Flt, Flt) = (v * 1e6).into();
fn loc(v: Meters) -> String {
let (x, y, z): (Flt, Flt, Flt) = (*v * 1e6).into();
format!("({}, {}, {}) um", x as i32, y as i32, z as i32)
}
/// M
pub struct Magnetization(pub Vec3);
pub struct Magnetization(pub Meters);
impl AbstractMeasurement for Magnetization {
fn eval(&self, state: &dyn GenericSim) -> String {
@@ -73,7 +75,7 @@ impl AbstractMeasurement for Magnetization {
}
/// B
pub struct MagneticFlux(pub Vec3);
pub struct MagneticFlux(pub Meters);
impl AbstractMeasurement for MagneticFlux {
fn eval(&self, state: &dyn GenericSim) -> String {
@@ -83,7 +85,7 @@ impl AbstractMeasurement for MagneticFlux {
}
/// H
pub struct MagneticStrength(pub Vec3);
pub struct MagneticStrength(pub Meters);
impl AbstractMeasurement for MagneticStrength {
fn eval(&self, state: &dyn GenericSim) -> String {
@@ -92,7 +94,7 @@ impl AbstractMeasurement for MagneticStrength {
}
}
pub struct ElectricField(pub Vec3);
pub struct ElectricField(pub Meters);
impl AbstractMeasurement for ElectricField {
fn eval(&self, state: &dyn GenericSim) -> String {

4
src/render.rs
View File

@@ -1,5 +1,5 @@
use ansi_term::Color::RGB;
use crate::geom::{Vec2, Vec3};
use crate::geom::{Meters, Vec2, Vec3};
use crate::{flt::{Flt, Real}, Material as _};
use crate::mat;
use crate::sim::{Cell, GenericSim};
@@ -92,7 +92,7 @@ impl<'a> RenderSteps<'a> {
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());
let z_m = self.z as Flt * self.sim.feature_size();
self.sim.sample(Vec3::new(x_m, y_m, z_m))
self.sim.sample(Meters(Vec3::new(x_m, y_m, z_m)))
}
////////////// Ex/Ey/Bz configuration ////////////

122
src/sim.rs
View File

@@ -1,5 +1,5 @@
use crate::{flt::{Flt, Real}, consts};
use crate::geom::{Coord, Vec3};
use crate::geom::{Coord, Index, Meters, Vec3, Vec3u};
use crate::mat::{self, GenericMaterial, Material};
use log::trace;
@@ -9,38 +9,61 @@ use std::convert::From;
use std::iter::Sum;
pub trait GenericSim {
fn sample(&self, pos_meters: Vec3) -> Cell<mat::Static>;
/// TODO: DEPRECATED
fn get(&self, at: Coord) -> Cell<mat::Static> {
let at_vec: Vec3 = at.into();
self.sample(at_vec * self.feature_size())
}
fn width(&self) -> u32;
fn height(&self) -> u32;
fn depth(&self) -> u32;
fn sample(&self, pos: Meters) -> Cell<mat::Static>;
fn impulse_e_meters(&mut self, pos: Meters, amount: Vec3);
fn size(&self) -> Index;
fn feature_size(&self) -> Flt;
fn time(&self) -> Flt;
fn timestep(&self) -> Flt;
fn step_no(&self) -> u64;
fn width(&self) -> u32 {
self.size().x()
}
fn height(&self) -> u32 {
self.size().y()
}
fn depth(&self) -> u32 {
self.size().z()
}
fn time(&self) -> Flt {
self.timestep() * self.step_no() as Flt
}
}
impl<'a> dyn GenericSim + 'a {
pub fn get<C: Coord>(&self, at: C) -> Cell<mat::Static> {
self.sample(at.to_meters(self.feature_size()))
}
/// Apply `F` to each Cell, and sum the results.
pub fn map_sum<F: Fn(&Cell<mat::Static>) -> R + , R: Sum<R>>(&self, f: F) -> R {
self.map_sum_enumerated(|_at, cell| f(cell))
self.map_sum_enumerated(|_at: Index, cell| f(cell))
}
pub fn map_sum_enumerated<F: Fn(Coord, &Cell<mat::Static>) -> R, R: Sum<R>>(&self, f: F) -> R {
pub fn map_sum_enumerated<C: Coord, F: Fn(C, &Cell<mat::Static>) -> R, R: Sum<R>>(&self, f: F) -> R {
let (w, h, d) = (self.width(), self.height(), self.depth());
(0..d).map(|z| (0..h).map(|y| (0..w).map(|x| {
let at = Coord::new(x, y, z);
f(at, &self.get(at))
let at = Index(Vec3u::new(x, y, z));
f(C::from_index(at, self.feature_size()), &self.get(at))
}).sum()).sum()).sum()
}
/// TODO: DEPRECATED
pub fn current(&self, c: Coord) -> Vec3 {
pub fn current<C: Coord>(&self, c: C) -> Vec3 {
self.get(c).current_density() * (self.feature_size() * self.feature_size())
}
pub fn impulse_e<C: Coord>(&mut self, pos: C, amt: Vec3) {
self.impulse_e_meters(pos.to_meters(self.feature_size()), amt)
}
pub fn impulse_ex<C: Coord>(&mut self, c: C, ex: Flt) {
self.impulse_e(c, Vec3::new(ex, 0.0, 0.0));
}
pub fn impulse_ey<C: Coord>(&mut self, c: C, ey: Flt) {
self.impulse_e(c, Vec3::new(0.0, ey, 0.0));
}
pub fn impulse_ez<C: Coord>(&mut self, c: C, ez: Flt) {
self.impulse_e(c, Vec3::new(0.0, 0.0, ez));
}
}
#[derive(Default)]
@@ -52,7 +75,7 @@ pub struct SimState<M=GenericMaterial> {
}
impl<M: Material + Default> SimState<M> {
pub fn new(size: Coord, feature_size: Flt) -> Self {
pub fn new(size: Index, feature_size: Flt) -> Self {
Self {
cells: Array3::default((size.z() as _, size.y() as _, size.x() as _)),
scratch: Array3::default((size.z() as _, size.y() as _, size.x() as _)),
@@ -100,23 +123,10 @@ impl<M: Material + Clone + Default + Send + Sync> SimState<M> {
}
}
impl<M: Material + Default + Send + Sync> SimState<M> {
/// Apply `F` to each Cell, and sum the results.
pub fn map_sum<F: Fn(&Cell<M>) -> R + Sync + Send, R: Sum<R> + Send>(&self, f: F) -> R {
self.cells.view().into_par_iter().map(f).sum()
}
pub fn map_sum_enumerated<F: Fn(Coord, &Cell<M>) -> R + Sync + Send, R: Sum<R> + Send>(&self, f: F) -> R {
Zip::from(ndarray::indices_of(&self.cells)).and(&self.cells).into_par_iter().map(|((z, y, x), c)| {
f(Coord::new(x as _, y as _, z as _), c)
}).sum()
}
}
impl<M: Material> GenericSim for SimState<M> {
fn sample(&self, pos_meters: Vec3) -> Cell<mat::Static> {
fn sample(&self, pos: Meters) -> Cell<mat::Static> {
// TODO: smarter sampling than nearest neighbor?
let pos_sim = pos_meters / self.feature_size();
let pos_sim = pos.to_index(self.feature_size());
let idx = [pos_sim.z() as usize, pos_sim.y() as _, pos_sim.x() as _];
match self.cells.get(idx) {
Some(cell) => Cell {
@@ -130,20 +140,22 @@ impl<M: Material> GenericSim for SimState<M> {
}
}
fn width(&self) -> u32 {
self.cells.shape()[2] as _
fn impulse_e_meters(&mut self, pos: Meters, amount: Vec3) {
self.get_mut(pos).state.e += amount;
}
fn height(&self) -> u32 {
self.cells.shape()[1] as _
}
fn depth(&self) -> u32 {
self.cells.shape()[0] as _
fn size(&self) -> Index {
Index(Vec3u::new(
self.cells.shape()[2] as _,
self.cells.shape()[1] as _,
self.cells.shape()[0] as _,
))
}
fn feature_size(&self) -> Flt {
self.feature_size.into()
}
fn time(&self) -> Flt {
self.timestep() * self.step_no as Flt
fn timestep(&self) -> Flt {
(self.feature_size / consts::real::C()).into()
}
fn step_no(&self) -> u64 {
self.step_no
@@ -151,31 +163,19 @@ impl<M: Material> GenericSim for SimState<M> {
}
impl<M: Material> SimState<M> {
pub fn impulse_bz(&mut self, c: Coord, bz: Flt) {
pub fn impulse_bz<C: Coord>(&mut self, c: C, bz: Flt) {
self.get_mut(c).impulse_bz(bz);
}
}
impl<M> SimState<M> {
pub fn impulse_ex(&mut self, c: Coord, ex: Flt) {
self.get_mut(c).state.e += Vec3::new(ex, 0.0, 0.0);
pub fn get<C: Coord>(&self, c: C) -> &Cell<M> {
let at = c.to_index(self.feature_size.into());
&self.cells[at.row_major_idx()]
}
pub fn impulse_ey(&mut self, c: Coord, ey: Flt) {
self.get_mut(c).state.e += Vec3::new(0.0, ey, 0.0);
}
pub fn impulse_ez(&mut self, c: Coord, ez: Flt) {
self.get_mut(c).state.e += Vec3::new(0.0, 0.0, ez);
}
pub fn get(&self, c: Coord) -> &Cell<M> {
&self.cells[c.row_major_idx()]
}
pub fn get_mut(&mut self, c: Coord) -> &mut Cell<M> {
&mut self.cells[c.row_major_idx()]
}
pub fn timestep(&self) -> Flt {
(self.feature_size / consts::real::C()).into()
pub fn get_mut<C: Coord>(&mut self, c: C) -> &mut Cell<M> {
let at = c.to_index(self.feature_size.into());
&mut self.cells[at.row_major_idx()]
}
}