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

Parameterize Vec3 over Real type. BROKEN TESTS

Browse Source 
sim::pml_boundary_conditions_{x,y,z} are broken, among others.
Unsure if this is new. I don't see an obvious cause yet.
This commit is contained in:
Colin
2021-06-07 14:23:07 -07:00
parent 84a92ad572
commit 9f6d92c5fc
16 changed files with 548 additions and 402 deletions
Download Patch File Download Diff File Expand all files Collapse all files

12
examples/minimal_torus.rs
View File

@@ -30,14 +30,14 @@ fn main() {
driver.set_steps_per_frame(400);
driver.set_steps_per_stim(1);
let base = "minimal_torus-6";
let ferro_region = Torus::new_xy(Meters((half_width, half_width, half_depth).into()), ferro_major, ferro_minor);
let drive_region = Torus::new_xz(Meters((half_width - ferro_major, half_width, half_depth).into()), wire_major, wire_minor);
let sense_region = Torus::new_xz(Meters((half_width + ferro_major, half_width, half_depth).into()), wire_major, wire_minor);
let ferro_region = Torus::new_xy(Meters::new(half_width, half_width, half_depth), ferro_major, ferro_minor);
let drive_region = Torus::new_xz(Meters::new(half_width - ferro_major, half_width, half_depth), wire_major, wire_minor);
let sense_region = Torus::new_xz(Meters::new(half_width + ferro_major, half_width, half_depth), wire_major, wire_minor);
let boundary_xy = half_width - ferro_major - ferro_minor - buffer;
let boundary_z = half_depth - wire_major - wire_minor - buffer;
let boundary_lower = Meters((boundary_xy, boundary_xy, boundary_z).into());
let boundary_upper = Meters((width - boundary_xy, width - boundary_xy, depth - boundary_z).into());
let boundary_lower = Meters::new(boundary_xy, boundary_xy, boundary_z);
let boundary_upper = Meters::new(width - boundary_xy, width - boundary_xy, depth - boundary_z);
println!("boundary: {}um; {}um", m_to_um(boundary_xy), m_to_um(boundary_z));
let boundary_region = InvertedRegion::new(Cube::new(boundary_lower, boundary_upper));
@@ -52,7 +52,7 @@ fn main() {
driver.fill_region(&drive_region, mat::db::conductor(conductivity).into());
driver.fill_region(&sense_region, mat::db::conductor(conductivity).into());
// driver.add_pml_boundary(Meters((boundary_xy, boundary_xy, boundary_z).into()));
driver.add_classical_boundary(Meters((boundary_xy, boundary_xy, boundary_z).into()));
driver.add_classical_boundary(Meters::new(boundary_xy, boundary_xy, boundary_z));
// J=\sigma E
// dJ/dt = \sigma dE/dT

18
examples/toroid25d.rs
View File

@@ -64,31 +64,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::MagnetizationAt(
Meters((half_width + ferro_inner_rad + 2.0*feat_size, half_width, half_depth).into())
Meters::new(half_width + ferro_inner_rad + 2.0*feat_size, half_width, half_depth)
));
driver.add_measurement(meas::MagneticFluxAt(
Meters((half_width + ferro_inner_rad + 2.0*feat_size, half_width, half_depth).into())
Meters::new(half_width + ferro_inner_rad + 2.0*feat_size, half_width, half_depth)
));
driver.add_measurement(meas::MagneticStrengthAt(
Meters((half_width + ferro_inner_rad + 2.0*feat_size, half_width, half_depth).into())
Meters::new(half_width + ferro_inner_rad + 2.0*feat_size, half_width, half_depth)
));
driver.add_measurement(meas::MagnetizationAt(
Meters((half_width + ferro_inner_rad + 1.0*feat_size, half_width, half_depth).into())
Meters::new(half_width + ferro_inner_rad + 1.0*feat_size, half_width, half_depth)
));
driver.add_measurement(meas::MagneticFluxAt(
Meters((half_width + ferro_inner_rad + 1.0*feat_size, half_width, half_depth).into())
Meters::new(half_width + ferro_inner_rad + 1.0*feat_size, half_width, half_depth)
));
driver.add_measurement(meas::MagneticStrengthAt(
Meters((half_width + ferro_inner_rad + 1.0*feat_size, half_width, half_depth).into())
Meters::new(half_width + ferro_inner_rad + 1.0*feat_size, half_width, half_depth)
));
driver.add_measurement(meas::MagnetizationAt(
Meters((half_width + 0.5 * (ferro_inner_rad + ferro_outer_rad), half_width, half_depth).into())
Meters::new(half_width + 0.5 * (ferro_inner_rad + ferro_outer_rad), half_width, half_depth)
));
driver.add_measurement(meas::MagneticFluxAt(
Meters((half_width + 0.5 * (ferro_inner_rad + ferro_outer_rad), half_width, half_depth).into())
Meters::new(half_width + 0.5 * (ferro_inner_rad + ferro_outer_rad), half_width, half_depth)
));
driver.add_measurement(meas::MagneticStrengthAt(
Meters((half_width + 0.5 * (ferro_inner_rad + ferro_outer_rad), half_width, half_depth).into())
Meters::new(half_width + 0.5 * (ferro_inner_rad + ferro_outer_rad), half_width, half_depth)
));
let center = Vec2::new(half_width as Flt, half_width as Flt);

29
examples/wrapped_torus.rs
View File

@@ -17,7 +17,8 @@ fn main() {
let peak_current = 5e6;
let current_duration = 0.1e-9; // half-wavelength of the sine wave
let current_break = 0.1e-9; // time between 'set' pulse and 'clear' pulse
let conductivity = 5.0e7;
let drive_conductivity = 5.0;
let sense_conductivity = 5.0;
let from_m = |m: Flt| (m/feat_size).round() as u32;
let m_to_um = |m: Flt| (m * 1e6).round() as u32;
@@ -34,28 +35,28 @@ fn main() {
let mut driver: Driver<mat::GenericMaterial> = Driver::new(size_px, feat_size);
driver.set_steps_per_frame(500);
// driver.set_steps_per_stim(10);
let base = "wrapped_torus-19";
let base = "wrapped_torus-21-high-input-resistance";
let ferro1_region = Torus::new_xy(Meters((q1_width, half_height, half_depth).into()), ferro_major, ferro_minor);
let drive1_region = Torus::new_xz(Meters((q1_width - ferro_major, half_height, half_depth).into()), wire_major, wire_minor);
let sense1_region = Torus::new_xz(Meters((q1_width + ferro_major, half_height, half_depth).into()), wire_major, wire_minor);
let ferro1_region = Torus::new_xy(Meters::new(q1_width, half_height, half_depth), ferro_major, ferro_minor);
let drive1_region = Torus::new_xz(Meters::new(q1_width - ferro_major, half_height, half_depth), wire_major, wire_minor);
let sense1_region = Torus::new_xz(Meters::new(q1_width + ferro_major, half_height, half_depth), wire_major, wire_minor);
driver.fill_region(&ferro1_region, mat::db::linear_iron().into());
driver.fill_region(&drive1_region, mat::db::conductor(conductivity).into());
driver.fill_region(&sense1_region, mat::db::conductor(conductivity).into());
driver.fill_region(&drive1_region, mat::db::conductor(drive_conductivity).into());
driver.fill_region(&sense1_region, mat::db::conductor(sense_conductivity).into());
let ferro2_region = Torus::new_xy(Meters((q3_width, half_height, half_depth).into()), ferro_major, ferro_minor);
let drive2_region = Torus::new_xz(Meters((q3_width - ferro_major, half_height, half_depth).into()), wire_major, wire_minor);
let sense2_region = Torus::new_xz(Meters((q3_width + ferro_major, half_height, half_depth).into()), wire_major, wire_minor);
let ferro2_region = Torus::new_xy(Meters::new(q3_width, half_height, half_depth), ferro_major, ferro_minor);
let drive2_region = Torus::new_xz(Meters::new(q3_width - ferro_major, half_height, half_depth), wire_major, wire_minor);
let sense2_region = Torus::new_xz(Meters::new(q3_width + ferro_major, half_height, half_depth), wire_major, wire_minor);
driver.fill_region(&ferro2_region, mat::db::ferroxcube_3r1().into());
driver.fill_region(&drive2_region, mat::db::conductor(conductivity).into());
driver.fill_region(&sense2_region, mat::db::conductor(conductivity).into());
driver.fill_region(&drive2_region, mat::db::conductor(drive_conductivity).into());
driver.fill_region(&sense2_region, mat::db::conductor(sense_conductivity).into());
let boundary_xy = q1_width - ferro_major - ferro_minor - buffer;
let boundary_z = half_depth - wire_major - wire_minor - buffer;
println!("boundary: {}um; {}um", m_to_um(boundary_xy), m_to_um(boundary_z));
driver.add_pml_boundary(Meters((boundary_xy, boundary_xy, boundary_z).into()));
driver.add_pml_boundary(Meters::new(boundary_xy, boundary_xy, boundary_z));
// J=\sigma E
// dJ/dt = \sigma dE/dT
@@ -63,7 +64,7 @@ fn main() {
// dE/dt = dI/dt / (A*\sigma)
// if I = k*sin(w t) then dE/dt = k*w sin(w t) / (A*\sigma)
// i.e. dE/dt is proportional to I/(A*\sigma), multiplied by w (or, divided by wavelength)
let peak_stim = peak_current/current_duration / (drive1_region.cross_section() * conductivity);
let peak_stim = peak_current/current_duration / (drive1_region.cross_section() * drive_conductivity);
let pos_wave = Sinusoid1::from_wavelength(peak_stim, current_duration * 2.0)
.half_cycle();

12
src/bin/pml_tuning.rs
View File

@@ -1,13 +1,15 @@
use coremem::*;
use coremem::geom::*;
use coremem::real::Real as _;
use coremem::stim::AbstractStimulus;
use rand::rngs::StdRng;
use rand::{Rng as _, SeedableRng as _};
fn energy<R: Region>(s: &dyn GenericSim, reg: &R) -> Flt {
0.5 * s.map_sum_over_enumerated(reg, |_pos: Index, cell| {
let e = Real::half() * s.map_sum_over_enumerated(reg, |_pos: Index, cell| {
cell.e().mag_sq()
})
});
Flt::from_primitive(e)
}
fn energy_now_and_then<R: Region>(state: &mut StaticSim, reg: &R, frames: u32) -> (Flt, Flt) {
@@ -31,7 +33,7 @@ impl<F: Fn(Index) -> (Vec3, Flt) + Sync> AbstractStimulus for PmlStim<F> {
let (e, hz) = (self.f)(pos.to_index(self.feat_size));
let sig_angle = angle*hz;
let sig = sig_angle.sin();
e*gate*sig
e*Real::from_primitive(gate*sig)
}
}
@@ -104,8 +106,8 @@ fn state_for_pml(size: Index, boundary: Index, feat_size: Flt, sc_coeff: Flt, co
let mut state = StaticSim::new(size, feat_size);
let timestep = state.timestep();
state.fill_boundary_using(boundary, |boundary_ness| {
let b = boundary_ness.elem_pow(pow);
let coord_stretch = b * (sc_coeff / timestep);
let b = boundary_ness.elem_pow(Real::from_primitive(pow));
let coord_stretch = b * Real::from_primitive(sc_coeff / timestep);
let conductivity = Vec3::unit() * (b.mag() * cond_coeff / timestep);
unimplemented!()
// Static {

15
src/driver.rs
View File

@@ -1,7 +1,8 @@
use crate::{flt::Flt, mat};
use crate::{flt, flt::Flt, mat};
use crate::geom::{Coord, Index, Meters, Region, Vec3, Vec3u};
use crate::mat::{GenericMaterial, Material, Pml};
use crate::meas::{self, AbstractMeasurement};
use crate::real::Real as _;
use crate::render::{self, MultiRenderer, Renderer};
use crate::sim::{GenericSim, SimState};
use crate::stim::AbstractStimulus;
@@ -193,8 +194,8 @@ impl<M: Material + From<mat::Pml>> Driver<M> {
pub fn add_pml_boundary<C: Coord>(&mut self, thickness: C) {
let timestep = self.state.timestep();
self.state.fill_boundary_using(thickness, |boundary_ness| {
let b = boundary_ness.elem_pow(3.0);
let conductivity = b * (0.5 / timestep);
let b = boundary_ness.elem_pow(flt::Real::three());
let conductivity = b * (flt::Real::half() / timestep);
Pml::new(conductivity).into()
});
}
@@ -204,11 +205,11 @@ impl<M: Material + From<mat::Conductor>> Driver<M> {
pub fn add_classical_boundary<C: Coord>(&mut self, thickness: C) {
let timestep = self.state.timestep();
self.state.fill_boundary_using(thickness, |boundary_ness| {
let b = boundary_ness.elem_pow(3.0);
let cond = b * (0.5 / timestep);
let b = boundary_ness.elem_pow(flt::Real::three());
let cond = b * (flt::Real::half() / timestep);
let iso_cond = cond.x() + cond.y() + cond.z();
let iso_conductor = mat::db::conductor(iso_cond);
let iso_conductor = mat::db::conductor(Flt::from_primitive(iso_cond));
iso_conductor.into()
});
}
@@ -223,7 +224,7 @@ struct StimuliAdapter {
impl AbstractStimulus for StimuliAdapter {
fn at(&self, t_sec: Flt, pos: Meters) -> Vec3 {
self.stim.at(t_sec, pos) * (self.frame_interval as Flt)
self.stim.at(t_sec, pos) * flt::Real::from_f64(self.frame_interval as f64)
}
}

155
src/geom/region.rs
View File

@@ -1,5 +1,6 @@
use crate::flt::{Flt, Real};
use crate::geom::{Meters, Vec2, Vec3};
use crate::real::Real as _;
use dyn_clone::{self, DynClone};
use serde::{Serialize, Deserialize};
@@ -98,7 +99,7 @@ impl Region for Torus {
// and they all give the same answer
Vec3::new(self.major_rad.into(), 0.0, 0.0)
} else {
p_on_plane.with_mag(self.major_rad.into_inner())
p_on_plane.with_mag(self.major_rad)
};
let distance_to_circle = (rel_p - q).mag();
distance_to_circle < self.minor_rad.into_inner()
@@ -140,22 +141,22 @@ impl Cube {
Self { lower, upper }
}
pub fn x_range(&self) -> Range<Flt> {
self.lower.x()..self.upper.x()
Flt::from_primitive(self.lower.x())..Flt::from_primitive(self.upper.x())
}
pub fn y_range(&self) -> Range<Flt> {
self.lower.y()..self.upper.y()
Flt::from_primitive(self.lower.y())..Flt::from_primitive(self.upper.y())
}
pub fn z_range(&self) -> Range<Flt> {
self.lower.z()..self.upper.z()
Flt::from_primitive(self.lower.z())..Flt::from_primitive(self.upper.z())
}
}
#[typetag::serde]
impl Region for Cube {
fn contains(&self, p: Meters) -> bool {
self.x_range().contains(&p.x()) &&
self.y_range().contains(&p.y()) &&
self.z_range().contains(&p.z())
self.x_range().contains(&Flt::from_primitive(p.x())) &&
self.y_range().contains(&Flt::from_primitive(p.y())) &&
self.z_range().contains(&Flt::from_primitive(p.z()))
}
}
@@ -214,103 +215,103 @@ mod test {
use super::*;
#[test]
fn torus_as_sphere() {
let tor = Torus::new_xy(Meters((0.0, 0.0, 0.0).into()), 0.0, 10.0);
assert!(tor.contains(Meters((0.0, 0.0, 0.0).into())));
assert!(tor.contains(Meters((5.0, 0.0, 0.0).into())));
let tor = Torus::new_xy(Meters::new(0.0, 0.0, 0.0), 0.0, 10.0);
assert!(tor.contains(Meters::new(0.0, 0.0, 0.0)));
assert!(tor.contains(Meters::new(5.0, 0.0, 0.0)));
assert!(tor.contains(Meters((9.0, 0.0, 0.0).into())));
assert!(tor.contains(Meters((-9.0, 0.0, 0.0).into())));
assert!(tor.contains(Meters((0.0, 9.0, 0.0).into())));
assert!(tor.contains(Meters((0.0, -9.0, 0.0).into())));
assert!(tor.contains(Meters((0.0, 0.0, 9.0).into())));
assert!(tor.contains(Meters((0.0, 0.0, -9.0).into())));
assert!(tor.contains(Meters::new(9.0, 0.0, 0.0)));
assert!(tor.contains(Meters::new(-9.0, 0.0, 0.0)));
assert!(tor.contains(Meters::new(0.0, 9.0, 0.0)));
assert!(tor.contains(Meters::new(0.0, -9.0, 0.0)));
assert!(tor.contains(Meters::new(0.0, 0.0, 9.0)));
assert!(tor.contains(Meters::new(0.0, 0.0, -9.0)));
assert!(tor.contains(Meters((5.0, 5.0, 5.0).into())));
assert!(tor.contains(Meters((-5.0, -5.0, -5.0).into())));
assert!(tor.contains(Meters::new(5.0, 5.0, 5.0)));
assert!(tor.contains(Meters::new(-5.0, -5.0, -5.0)));
// boundary
assert!(!tor.contains(Meters((10.0, 0.0, 0.0).into())));
assert!(!tor.contains(Meters((-10.0, 0.0, 0.0).into())));
assert!(!tor.contains(Meters((0.0, 10.0, 0.0).into())));
assert!(!tor.contains(Meters((0.0, -10.0, 0.0).into())));
assert!(!tor.contains(Meters((0.0, 0.0, 10.0).into())));
assert!(!tor.contains(Meters((0.0, 0.0, -10.0).into())));
assert!(!tor.contains(Meters::new(10.0, 0.0, 0.0)));
assert!(!tor.contains(Meters::new(-10.0, 0.0, 0.0)));
assert!(!tor.contains(Meters::new(0.0, 10.0, 0.0)));
assert!(!tor.contains(Meters::new(0.0, -10.0, 0.0)));
assert!(!tor.contains(Meters::new(0.0, 0.0, 10.0)));
assert!(!tor.contains(Meters::new(0.0, 0.0, -10.0)));
assert!(!tor.contains(Meters((11.0, 0.0, 0.0).into())));
assert!(!tor.contains(Meters((9.0, 9.0, 9.0).into())));
assert!(!tor.contains(Meters((-9.0, -9.0, -9.0).into())));
assert!(!tor.contains(Meters((-9.0, 9.0, 0.0).into())));
assert!(!tor.contains(Meters::new(11.0, 0.0, 0.0)));
assert!(!tor.contains(Meters::new(9.0, 9.0, 9.0)));
assert!(!tor.contains(Meters::new(-9.0, -9.0, -9.0)));
assert!(!tor.contains(Meters::new(-9.0, 9.0, 0.0)));
}
#[test]
fn torus_xy() {
let tor = Torus::new_xy(Meters((0.0, 0.0, 0.0).into()), 10.0, 2.0);
assert!(!tor.contains(Meters((0.0, 0.0, 0.0).into())));
let tor = Torus::new_xy(Meters::new(0.0, 0.0, 0.0), 10.0, 2.0);
assert!(!tor.contains(Meters::new(0.0, 0.0, 0.0)));
assert!(tor.contains(Meters((10.0, 0.0, 0.0).into())));
assert!(tor.contains(Meters((-10.0, 0.0, 0.0).into())));
assert!(tor.contains(Meters((0.0, 10.0, 0.0).into())));
assert!(tor.contains(Meters((0.0, -10.0, 0.0).into())));
assert!(!tor.contains(Meters((0.0, 0.0, 10.0).into())));
assert!(!tor.contains(Meters((0.0, 0.0, -10.0).into())));
assert!(tor.contains(Meters::new(10.0, 0.0, 0.0)));
assert!(tor.contains(Meters::new(-10.0, 0.0, 0.0)));
assert!(tor.contains(Meters::new(0.0, 10.0, 0.0)));
assert!(tor.contains(Meters::new(0.0, -10.0, 0.0)));
assert!(!tor.contains(Meters::new(0.0, 0.0, 10.0)));
assert!(!tor.contains(Meters::new(0.0, 0.0, -10.0)));
assert!(tor.contains(Meters((11.0, 0.0, 0.0).into())));
assert!(tor.contains(Meters((-11.0, 0.0, 0.0).into())));
assert!(tor.contains(Meters((11.0, 2.0, 0.0).into())));
assert!(tor.contains(Meters((11.0, 0.0, 1.0).into())));
assert!(tor.contains(Meters((-11.0, 0.0, 1.0).into())));
assert!(tor.contains(Meters((-7.0, 7.0, 1.0).into())));
assert!(tor.contains(Meters((-7.0, -7.0, -1.0).into())));
assert!(tor.contains(Meters::new(11.0, 0.0, 0.0)));
assert!(tor.contains(Meters::new(-11.0, 0.0, 0.0)));
assert!(tor.contains(Meters::new(11.0, 2.0, 0.0)));
assert!(tor.contains(Meters::new(11.0, 0.0, 1.0)));
assert!(tor.contains(Meters::new(-11.0, 0.0, 1.0)));
assert!(tor.contains(Meters::new(-7.0, 7.0, 1.0)));
assert!(tor.contains(Meters::new(-7.0, -7.0, -1.0)));
// boundary
assert!(!tor.contains(Meters((12.0, 0.0, 0.0).into())));
assert!(!tor.contains(Meters((10.0, 0.0, 2.0).into())));
assert!(!tor.contains(Meters((-10.0, 0.0, 2.0).into())));
assert!(!tor.contains(Meters::new(12.0, 0.0, 0.0)));
assert!(!tor.contains(Meters::new(10.0, 0.0, 2.0)));
assert!(!tor.contains(Meters::new(-10.0, 0.0, 2.0)));
}
#[test]
fn torus_xz() {
let tor = Torus::new_xz(Meters((0.0, 0.0, 0.0).into()), 10.0, 2.0);
assert!(!tor.contains(Meters((0.0, 0.0, 0.0).into())));
let tor = Torus::new_xz(Meters::new(0.0, 0.0, 0.0), 10.0, 2.0);
assert!(!tor.contains(Meters::new(0.0, 0.0, 0.0)));
assert!(tor.contains(Meters((10.0, 0.0, 0.0).into())));
assert!(tor.contains(Meters((-10.0, 0.0, 0.0).into())));
assert!(!tor.contains(Meters((0.0, 10.0, 0.0).into())));
assert!(!tor.contains(Meters((0.0, -10.0, 0.0).into())));
assert!(tor.contains(Meters((0.0, 0.0, 10.0).into())));
assert!(tor.contains(Meters((0.0, 0.0, -10.0).into())));
assert!(tor.contains(Meters::new(10.0, 0.0, 0.0)));
assert!(tor.contains(Meters::new(-10.0, 0.0, 0.0)));
assert!(!tor.contains(Meters::new(0.0, 10.0, 0.0)));
assert!(!tor.contains(Meters::new(0.0, -10.0, 0.0)));
assert!(tor.contains(Meters::new(0.0, 0.0, 10.0)));
assert!(tor.contains(Meters::new(0.0, 0.0, -10.0)));
assert!(tor.contains(Meters((11.0, 0.0, 0.0).into())));
assert!(tor.contains(Meters((-11.0, 0.0, 0.0).into())));
assert!(tor.contains(Meters((11.0, 1.0, 0.0).into())));
assert!(tor.contains(Meters((11.0, 0.0, 2.0).into())));
assert!(tor.contains(Meters((-11.0, 1.0, 0.0).into())));
assert!(tor.contains(Meters((-7.0, 1.0, 7.0).into())));
assert!(tor.contains(Meters((-7.0, -1.0, -7.0).into())));
assert!(tor.contains(Meters::new(11.0, 0.0, 0.0)));
assert!(tor.contains(Meters::new(-11.0, 0.0, 0.0)));
assert!(tor.contains(Meters::new(11.0, 1.0, 0.0)));
assert!(tor.contains(Meters::new(11.0, 0.0, 2.0)));
assert!(tor.contains(Meters::new(-11.0, 1.0, 0.0)));
assert!(tor.contains(Meters::new(-7.0, 1.0, 7.0)));
assert!(tor.contains(Meters::new(-7.0, -1.0, -7.0)));
// boundary
assert!(!tor.contains(Meters((12.0, 0.0, 0.0).into())));
assert!(!tor.contains(Meters((10.0, 2.0, 0.0).into())));
assert!(!tor.contains(Meters((-10.0, 2.0, 0.0).into())));
assert!(!tor.contains(Meters::new(12.0, 0.0, 0.0)));
assert!(!tor.contains(Meters::new(10.0, 2.0, 0.0)));
assert!(!tor.contains(Meters::new(-10.0, 2.0, 0.0)));
}
#[test]
fn torus_with_center() {
let tor = Torus::new_xy(Meters((1.0, 2.0, 3.0).into()), 0.0, 10.0);
assert!(tor.contains(Meters((1.0, 2.0, 3.0).into())));
assert!(tor.contains(Meters((1.0, 2.0, 12.0).into())));
assert!(!tor.contains(Meters((1.0, 2.0, 13.0).into())));
let tor = Torus::new_xy(Meters::new(1.0, 2.0, 3.0), 0.0, 10.0);
assert!(tor.contains(Meters::new(1.0, 2.0, 3.0)));
assert!(tor.contains(Meters::new(1.0, 2.0, 12.0)));
assert!(!tor.contains(Meters::new(1.0, 2.0, 13.0)));
let tor = Torus::new_xy(Meters((1.0, 2.0, 3.0).into()), 10.0, 3.0);
assert!(!tor.contains(Meters((1.0, 2.0, 3.0).into())));
assert!(tor.contains(Meters((1.0, 14.0, 3.0).into())));
assert!(!tor.contains(Meters((1.0, 15.0, 3.0).into())));
let tor = Torus::new_xy(Meters::new(1.0, 2.0, 3.0), 10.0, 3.0);
assert!(!tor.contains(Meters::new(1.0, 2.0, 3.0)));
assert!(tor.contains(Meters::new(1.0, 14.0, 3.0)));
assert!(!tor.contains(Meters::new(1.0, 15.0, 3.0)));
let tor = Torus::new_xz(Meters((1.0, 2.0, 3.0).into()), 5.0, 5.0);
assert!(!tor.contains(Meters((1.0, 2.0, 3.0).into())));
assert!(tor.contains(Meters((2.0, 2.0, 3.0).into())));
assert!(tor.contains(Meters((1.0, 2.0, 12.0).into())));
assert!(!tor.contains(Meters((1.0, 2.0, 13.0).into())));
let tor = Torus::new_xz(Meters::new(1.0, 2.0, 3.0), 5.0, 5.0);
assert!(!tor.contains(Meters::new(1.0, 2.0, 3.0)));
assert!(tor.contains(Meters::new(2.0, 2.0, 3.0)));
assert!(tor.contains(Meters::new(1.0, 2.0, 12.0)));
assert!(!tor.contains(Meters::new(1.0, 2.0, 13.0)));
}
}

13
src/geom/units.rs
View File

@@ -1,4 +1,5 @@
use crate::flt::Flt;
use crate::flt::{Flt, Real};
use crate::real::{Real as _, ToFloat};
use serde::{Serialize, Deserialize};
use super::{Vec3, Vec3u};
use std::fmt::{self, Display};
@@ -15,12 +16,18 @@ pub trait Coord: Copy + Clone {
#[derive(Copy, Clone, Debug, Serialize, Deserialize)]
pub struct Meters(pub Vec3);
impl Meters {
pub fn new<R: ToFloat>(x: R, y: R, z: R) -> Self {
Self(Vec3::new(x, y, z).into())
}
}
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).round().into())
Index((self.0 / Real::from_f64(feature_size)).round().into())
}
fn from_meters(other: Meters, _feature_size: Flt) -> Self {
other
@@ -74,7 +81,7 @@ impl Index {
impl Coord for Index {
fn to_meters(&self, feature_size: Flt) -> Meters {
Meters(Vec3::from(self.0) * feature_size)
Meters(Vec3::from(self.0) * Real::from_f64(feature_size))
}
fn to_index(&self, _feature_size: Flt) -> Index {
*self

279
src/geom/vec.rs
View File

@@ -1,5 +1,5 @@
use crate::flt::{self, Flt};
use crate::real::Real;
use crate::real::{Real, ToFloat};
use super::Vec3u;
use serde::{Serialize, Deserialize};
@@ -63,11 +63,11 @@ impl<R> Into<(R, R)> for Vec2<R> {
}
}
impl<R> Vec2<R> {
pub fn new<R2: Into<R>>(x: R2, y: R2) -> Self {
impl<R: Real> Vec2<R> {
pub fn new<R2: ToFloat>(x: R2, y: R2) -> Self {
Self {
x: x.into(),
y: y.into(),
x: R::from_primitive(x),
y: R::from_primitive(y),
}
}
}
@@ -124,71 +124,113 @@ impl<R: Real> Vec2<R> {
}
#[derive(Copy, Clone, Debug, Default, PartialEq, Serialize, Deserialize)]
pub struct Vec3 {
pub(crate) x: flt::Real,
pub(crate) y: flt::Real,
pub(crate) z: flt::Real,
pub struct Vec3<R=flt::Real> {
pub(crate) x: R,
pub(crate) y: R,
pub(crate) z: R,
}
impl Vec3 {
pub fn new(x: Flt, y: Flt, z: Flt) -> Self {
Self {
x: x.into(),
y: y.into(),
z: z.into()
}
}
pub fn to_tuple(self) -> (Flt, Flt, Flt) {
impl<R> Vec3<R> {
pub fn to_tuple(self) -> (R, R, R) {
self.into()
}
pub fn unit() -> Self {
Self::new(1.0, 1.0, 1.0)
}
pub fn unit_x() -> Self {
Self::new(1.0, 0.0, 0.0)
}
pub fn unit_y() -> Self {
Self::new(0.0, 1.0, 0.0)
}
pub fn unit_z() -> Self {
Self::new(0.0, 0.0, 1.0)
}
pub fn zero() -> Self {
Self::new(0.0, 0.0, 0.0)
}
pub fn x(&self) -> Flt {
self.x.into()
}
pub fn y(&self) -> Flt {
self.y.into()
}
pub fn z(&self) -> Flt {
self.z.into()
}
pub fn set_x(&mut self, x: Flt) {
self.x = flt::Real::from_inner(x);
}
pub fn set_y(&mut self, y: Flt) {
self.y = flt::Real::from_inner(y);
}
pub fn set_z(&mut self, z: Flt) {
self.z = flt::Real::from_inner(z);
}
pub fn xy(&self) -> Vec2 {
Vec2::new(self.x(), self.y())
}
pub fn mag(&self) -> Flt {
self.mag_sq().sqrt()
}
pub fn mag_sq(&self) -> Flt {
let Vec3 { x, y, z } = *self;
(x*x + y*y + z*z).into_inner()
}
pub fn component_sum(&self) -> Flt {
(self.x + self.y + self.z).into()
}
pub fn dot(&self, other: Self) -> Flt {
impl<R: Real> Vec3<R> {
pub fn new<R2: ToFloat>(x: R2, y: R2, z: R2) -> Self {
Self {
x: R::from_primitive(x),
y: R::from_primitive(y),
z: R::from_primitive(z),
}
}
pub fn new_x<R2: ToFloat>(x: R2) -> Self {
Self {
x: R::from_primitive(x),
y: R::zero(),
z: R::zero(),
}
}
pub fn new_y<R2: ToFloat>(y: R2) -> Self {
Self {
x: R::zero(),
y: R::from_primitive(y),
z: R::zero(),
}
}
pub fn new_z<R2: ToFloat>(z: R2) -> Self {
Self {
x: R::zero(),
y: R::zero(),
z: R::from_primitive(z),
}
}
pub fn uniform<R2: ToFloat>(u: R2) -> Self {
let u = R::from_primitive(u);
Self {
x: u,
y: u,
z: u,
}
}
pub fn cast<R2: Real>(&self) -> Vec3<R2> {
Vec3::new(
R2::from_primitive(self.x),
R2::from_primitive(self.y),
R2::from_primitive(self.z)
)
}
pub fn unit() -> Self {
Self::new(R::one(), R::one(), R::one())
}
pub fn unit_x() -> Self {
Self::new(R::one(), R::zero(), R::zero())
}
pub fn unit_y() -> Self {
Self::new(R::zero(), R::one(), R::zero())
}
pub fn unit_z() -> Self {
Self::new(R::zero(), R::zero(), R::one())
}
pub fn zero() -> Self {
Self::new(R::zero(), R::zero(), R::zero())
}
pub fn set_x<R2: ToFloat>(&mut self, x: R2) {
self.x = R::from_primitive(x);
}
pub fn set_y<R2: ToFloat>(&mut self, y: R2) {
self.y = R::from_primitive(y);
}
pub fn set_z<R2: ToFloat>(&mut self, z: R2) {
self.z = R::from_primitive(z);
}
pub fn x(&self) -> R {
self.x
}
pub fn y(&self) -> R {
self.y
}
pub fn z(&self) -> R {
self.z
}
pub fn xy(&self) -> Vec2<R> {
Vec2::new(self.x(), self.y())
}
pub fn mag(&self) -> R {
self.mag_sq().sqrt()
}
pub fn mag_sq(&self) -> R {
let Vec3 { x, y, z } = *self;
x*x + y*y + z*z
}
pub fn component_sum(&self) -> R {
self.x + self.y + self.z
}
pub fn dot(&self, other: Self) -> R {
self.elem_mul(other).component_sum()
}
@@ -211,7 +253,7 @@ impl Vec3 {
}
}
pub fn elem_pow(&self, pow: Flt) -> Self {
pub fn elem_pow(&self, pow: R) -> Self {
Self::new(
self.x().powf(pow),
self.y().powf(pow),
@@ -230,61 +272,60 @@ impl Vec3 {
/// Return (x^-1, y^-1, z^-1)
pub fn elem_inv(&self) -> Self {
Self::new(1.0 / self.x(), 1.0 / self.y(), 1.0 / self.z())
let one = R::one();
Self::new(one / self.x(), one / self.y(), one / self.z())
}
pub fn exp(&self) -> Self {
Self::new(self.x().exp(), self.y().exp(), self.z().exp())
}
pub fn with_mag(&self, new_mag: Flt) -> Vec3 {
if new_mag == 0.0 {
pub fn with_mag(&self, new_mag: R) -> Self {
if new_mag.is_zero() {
// avoid div-by-zero if self.mag() == 0 and new_mag == 0
Self::zero()
} else {
let scale = flt::Real::from_inner(new_mag) / self.mag();
self.clone() * scale.into_inner()
let scale = new_mag / self.mag();
*self * scale
}
}
/// Normalizes to 1.0 length, unless it's the zero vector, in which case length remains zero.
pub fn norm(&self) -> Vec3 {
pub fn norm(&self) -> Self {
if *self == Self::zero() {
*self
} else {
self.with_mag(1.0)
self.with_mag(R::one())
}
}
pub fn round(&self) -> Vec3 {
pub fn round(&self) -> Self {
Self::new(self.x().round(), self.y().round(), self.z().round())
}
}
impl Into<(Flt, Flt, Flt)> for Vec3 {
fn into(self) -> (Flt, Flt, Flt) {
(self.x(), self.y(), self.z())
impl<R> Into<(R, R, R)> for Vec3<R> {
fn into(self) -> (R, R, R) {
(self.x, self.y, self.z)
}
}
impl From<(Flt, Flt, Flt)> for Vec3 {
fn from((x, y, z): (Flt, Flt, Flt)) -> Self {
Self::new(x, y, z)
}
}
impl From<(flt::Real, flt::Real, flt::Real)> for Vec3 {
fn from((x, y, z): (flt::Real, flt::Real, flt::Real)) -> Self {
impl<R> From<(R, R, R)> for Vec3<R> {
fn from((x, y, z): (R, R, R)) -> Self {
Self { x, y, z }
}
}
impl From<Vec3u> for Vec3 {
impl<R: Real> From<Vec3u> for Vec3<R> {
fn from(v: Vec3u) -> Self {
Self::new(v.x() as _, v.y() as _, v.z() as _)
Self::new(
R::from_primitive(v.x()),
R::from_primitive(v.y()),
R::from_primitive(v.z())
)
}
}
impl AddAssign for Vec3 {
impl<R: Real> AddAssign for Vec3<R> {
fn add_assign(&mut self, other: Self) {
self.x += other.x;
self.y += other.y;
@@ -292,7 +333,7 @@ impl AddAssign for Vec3 {
}
}
impl Add for Vec3 {
impl<R: Real> Add for Vec3<R> {
type Output = Self;
fn add(self, other: Self) -> Self {
let mut work = self.clone();
@@ -301,7 +342,7 @@ impl Add for Vec3 {
}
}
impl Sub for Vec3 {
impl<R: Real> Sub for Vec3<R> {
type Output = Self;
fn sub(self, other: Self) -> Self {
Self {
@@ -312,7 +353,7 @@ impl Sub for Vec3 {
}
}
impl Neg for Vec3 {
impl<R: Real> Neg for Vec3<R> {
type Output = Self;
fn neg(self) -> Self {
Self {
@@ -323,77 +364,45 @@ impl Neg for Vec3 {
}
}
impl MulAssign<flt::Real> for Vec3 {
fn mul_assign(&mut self, other: flt::Real) {
impl<R: Real> MulAssign<R> for Vec3<R> {
fn mul_assign(&mut self, other: R) {
self.x *= other;
self.y *= other;
self.z *= other;
}
}
impl Mul<flt::Real> for Vec3 {
impl<R: Real> Mul<R> for Vec3<R> {
type Output = Self;
fn mul(self, other: flt::Real) -> Self {
fn mul(self, other: R) -> Self {
let mut work = self.clone();
work *= other;
work
}
}
impl MulAssign<Flt> for Vec3 {
fn mul_assign(&mut self, other: Flt) {
self.x *= other;
self.y *= other;
self.z *= other;
impl<R: Real> DivAssign<R> for Vec3<R> {
fn div_assign(&mut self, other: R) {
*self *= R::one() / other;
}
}
impl Mul<Flt> for Vec3 {
impl<R: Real> Div<R> for Vec3<R> {
type Output = Self;
fn mul(self, other: Flt) -> Self {
let mut work = self.clone();
work *= other;
work
}
}
impl DivAssign<flt::Real> for Vec3 {
fn div_assign(&mut self, other: flt::Real) {
*self *= flt::Real::from_inner(1.0) / other;
}
}
impl Div<flt::Real> for Vec3 {
type Output = Self;
fn div(self, other: flt::Real) -> Self {
fn div(self, other: R) -> Self {
let mut work = self.clone();
work /= other;
work
}
}
impl DivAssign<Flt> for Vec3 {
fn div_assign(&mut self, other: Flt) {
*self *= 1.0 / other;
}
}
impl Div<Flt> for Vec3 {
type Output = Self;
fn div(self, other: Flt) -> Self {
let mut work = self.clone();
work /= other;
work
}
}
impl Sum for Vec3 {
impl<R: Real> Sum for Vec3<R> {
fn sum<I: Iterator<Item = Self>>(iter: I) -> Self {
iter.fold(Self::default(), |a, b| a + b)
iter.fold(Self::zero(), |a, b| a + b)
}
}
impl fmt::LowerExp for Vec3 {
impl<R: Real> fmt::LowerExp for Vec3<R> {
fn fmt(&self, f: &mut fmt::Formatter<'_>) -> Result<(), fmt::Error> {
write!(f, "(")?;
fmt::LowerExp::fmt(&self.x(), f)?;
@@ -405,7 +414,7 @@ impl fmt::LowerExp for Vec3 {
}
}
impl fmt::Display for Vec3 {
impl<R: Real> fmt::Display for Vec3<R> {
fn fmt(&self, f: &mut fmt::Formatter<'_>) -> Result<(), fmt::Error> {
write!(f, "(")?;
fmt::Display::fmt(&self.x(), f)?;

5
src/geom/vecu.rs
View File

@@ -1,4 +1,5 @@
use super::Vec3;
use crate::geom::Vec3;
use crate::real::{Real as _, ToFloat as _};
use serde::{Serialize, Deserialize};
use std::fmt::{self, Display};
@@ -41,7 +42,7 @@ impl From<(u32, u32, u32)> for Vec3u {
impl From<Vec3> for Vec3u {
fn from(v: Vec3) -> Self {
Self::new(v.x() as _, v.y() as _, v.z() as _)
Self::new(v.x().to_f64() as _, v.y().to_f64() as _, v.z().to_f64() as _)
}
}

9
src/lib.rs
View File

@@ -33,8 +33,8 @@ pub(crate) mod flt {
use decorum::{R32, R64};
// Total enforces zero constraints (allows NaN, INFs).
// decorum::Real forbids NaN and INFs
pub(crate) type T64 = decorum::Total<f64>;
pub(crate) type T32 = decorum::Total<f32>;
pub type T64 = decorum::Total<f64>;
pub type T32 = decorum::Total<f32>;
// Use checked floats for debug and fast floats for release
#[cfg(debug)]
@@ -42,7 +42,7 @@ pub(crate) mod flt {
use super::*;
//pub(crate) type Real = R32;
//pub type Flt = f32;
pub(crate) type Real = R64;
pub type Real = R64;
pub type Flt = f64;
}
#[cfg(not(debug))]
@@ -51,12 +51,13 @@ pub(crate) mod flt {
// pub(crate) type Real = T32;
// pub type Flt = f32;
// pub(crate) type Real = T64;
pub(crate) type Real = R64;
pub type Real = R64;
pub type Flt = f64;
}
pub use defaults::*;
}
pub use flt::Flt;
pub use flt::Real;
#[allow(non_snake_case, unused)]
pub mod consts {

23
src/mat.rs
View File

@@ -1,6 +1,7 @@
use crate::{CellState, consts};
use crate::flt::{Flt, Real};
use crate::geom::{Line2d, Vec2, Vec3, Polygon2d};
use crate::real::Real as _;
use crate::sim::{PmlParameters, PmlState, StepParameters, StepParametersMut};
use lazy_static::lazy_static;
@@ -95,7 +96,7 @@ pub struct Conductor {
impl Conductor {
pub fn new(conductivity: Flt) -> Self {
Self::new_anisotropic(Vec3::unit() * conductivity)
Self::new_anisotropic(Vec3::unit() * Real::from_f64(conductivity))
}
pub fn new_anisotropic(conductivity: Vec3) -> Self {
Self {
@@ -120,7 +121,7 @@ pub struct LinearMagnet {
impl LinearMagnet {
pub fn new(relative_permeability: Flt) -> Self {
Self::new_anisotropic(Vec3::unit() * relative_permeability)
Self::new_anisotropic(Vec3::unit() * Real::from_f64(relative_permeability))
}
pub fn new_anisotropic(relative_permeability: Vec3) -> Self {
Self { relative_permeability, m: Vec3::zero() }
@@ -182,9 +183,21 @@ impl<T: PiecewiseLinearFerromagnet> Material for T {
let target_hm = h + m + delta_b * consts::real::MU0_INV();
// TODO: this is probably not the best way to generalize a BH curve into 3d.
let (_hx, mx) = mh_curve.move_to(h.x(), m.x(), target_hm.x());
let (_hy, my) = mh_curve.move_to(h.y(), m.y(), target_hm.y());
let (_hz, mz) = mh_curve.move_to(h.z(), m.z(), target_hm.z());
let (_hx, mx) = mh_curve.move_to(
Flt::from_primitive(h.x()),
Flt::from_primitive(m.x()),
Flt::from_primitive(target_hm.x())
);
let (_hy, my) = mh_curve.move_to(
Flt::from_primitive(h.y()),
Flt::from_primitive(m.y()),
Flt::from_primitive(target_hm.y())
);
let (_hz, mz) = mh_curve.move_to(
Flt::from_primitive(h.z()),
Flt::from_primitive(m.z()),
Flt::from_primitive(target_hm.z())
);
*self.m_mut() = Vec3::new(mx, my, mz);
// let ret = Vec3::new(hx, hy, hz);

37
src/meas.rs
View File

@@ -1,6 +1,7 @@
use crate::flt::Flt;
use crate::flt::{Flt, Real};
use crate::geom::{Meters, Region, Torus, Vec3, WorldRegion};
use crate::mat::Material as _;
use crate::real::Real as _;
use crate::sim::GenericSim;
use common_macros::b_tree_map;
use dyn_clone::{self, DynClone};
@@ -84,19 +85,19 @@ impl Current {
region: Box::new(r)
}
}
fn data(&self, state: &dyn GenericSim) -> (Flt, Vec3) {
fn data(&self, state: &dyn GenericSim) -> (Real, Vec3) {
let FieldSample(volume, current_mag, current_vec) = state.map_sum_over_enumerated(&*self.region, |coord: Meters, _cell| {
let current = state.current(coord);
FieldSample(1, current.mag(), current)
});
let mean_current_mag = current_mag / (volume as Flt);
let mean_current_vec = current_vec / (volume as Flt);
let mean_current_mag = current_mag / (Real::from_u32(volume));
let mean_current_vec = current_vec / (Real::from_u32(volume));
(mean_current_mag, mean_current_vec)
}
}
#[derive(Default)]
struct FieldSample(usize, Flt, Vec3);
struct FieldSample(u32, Real, Vec3);
impl std::iter::Sum for FieldSample {
fn sum<I>(iter: I) -> Self
@@ -172,7 +173,7 @@ impl CurrentLoop {
region: r,
}
}
fn data(&self, state: &dyn GenericSim) -> Flt {
fn data(&self, state: &dyn GenericSim) -> Real {
let FieldSample(volume, directed_current, _current_vec) = state.map_sum_over_enumerated(&self.region, |coord: Meters, _cell| {
let normal = self.region.axis();
let to_coord = *coord - *self.region.center();
@@ -181,7 +182,7 @@ impl CurrentLoop {
let directed_current = current.dot(tangent);
FieldSample(1, directed_current, current)
});
let mean_directed_current = directed_current / (volume as Flt);
let mean_directed_current = directed_current / Real::from_u32(volume);
let cross_section = self.region.cross_section() / (state.feature_size() * state.feature_size());
let cross_sectional_current = mean_directed_current * cross_section;
cross_sectional_current
@@ -217,7 +218,7 @@ impl MagneticLoop {
region: r,
}
}
fn data(&self, state: &dyn GenericSim) -> (Flt, Flt, Flt) {
fn data(&self, state: &dyn GenericSim) -> (Real, Real, Real) {
let FieldSamples([
FieldSample(volume, directed_m, _m_vec),
FieldSample(_, directed_b, _b_vec),
@@ -241,11 +242,11 @@ impl MagneticLoop {
])
});
// let cross_section = self.region.cross_section() / (state.feature_size() * state.feature_size());
let mean_directed_m = directed_m / (volume as Flt);
let mean_directed_m = directed_m / Real::from_u32(volume);
// let cross_sectional_m = mean_directed_m * cross_section;
let mean_directed_b = directed_b / (volume as Flt);
let mean_directed_b = directed_b / Real::from_u32(volume);
// let cross_sectional_b = mean_directed_b * cross_section;
let mean_directed_h = directed_h / (volume as Flt);
let mean_directed_h = directed_h / Real::from_u32(volume);
// let cross_sectional_h = mean_directed_h * cross_section;
// format!(
// "M({}): {:.2e}; B({}): {:.2e}; H({}): {:.2e}",
@@ -298,7 +299,7 @@ impl MagneticFlux {
let mag = b.mag();
FieldSample(1, mag, b)
});
let mean_mag = mag_vec / (volume as Flt);
let mean_mag = mag_vec / Real::from_u32(volume);
mean_mag
}
}
@@ -337,7 +338,7 @@ impl Magnetization {
let mag = m.mag();
FieldSample(1, mag, m)
});
let mean_mag = mag_vec / (volume as Flt);
let mean_mag = mag_vec / Real::from_u32(volume);
mean_mag
}
}
@@ -357,7 +358,7 @@ impl AbstractMeasurement for Magnetization {
}
fn loc(v: Meters) -> String {
format!("{:.0} um", *v * 1e6)
format!("{:.0} um", *v * Real::from_u64(1_000_000))
}
/// M
@@ -447,7 +448,7 @@ impl Energy {
region: Box::new(region),
}
}
fn data(&self, state: &dyn GenericSim) -> Flt {
fn data(&self, state: &dyn GenericSim) -> Real {
// Potential energy stored in a E/M field:
// https://en.wikipedia.org/wiki/Magnetic_energy
// https://en.wikipedia.org/wiki/Electric_potential_energy#Energy_stored_in_an_electrostatic_field_distribution
@@ -456,7 +457,7 @@ impl Energy {
// U(E) = 1/2 \int E . D dV
#[allow(non_snake_case)]
let dV = state.feature_volume();
let e = 0.5 * dV * state.map_sum_over(&*self.region, |cell| {
let e = Real::from_primitive(0.5 * dV) * state.map_sum_over(&*self.region, |cell| {
// E . D = E . (E + P) = E.E since we don't model polarization fields
cell.h().dot(cell.b()) + cell.e().mag_sq()
});
@@ -498,11 +499,11 @@ impl Power {
// Power is P = IV = A*J*V = L^2*J.(LE) = L^3 J.E
// where L is feature size.
#[allow(non_snake_case)]
let dV = state.feature_volume();
let dV = Real::from_primitive(state.feature_volume());
let power = dV * state.map_sum_over(&*self.region, |cell| {
cell.current_density().dot(cell.e())
});
power
Flt::from_primitive(power)
}
}

130
src/real.rs
View File

@@ -1,9 +1,10 @@
use std::ops::AddAssign;
use std::fmt;
use std::ops::{AddAssign, MulAssign};
use decorum::cmp::IntrinsicOrd;
pub use decorum::Real as decorum_Real;
pub use num::Zero;
/// This exists to allow configuration over # of bits (f32 v.s. f64) as well as
/// constraints.
pub trait Real: decorum::Real + IntrinsicOrd + AddAssign {
pub trait ToFloat {
fn to_f32(&self) -> f32 {
self.to_f64() as _
}
@@ -12,39 +13,142 @@ pub trait Real: decorum::Real + IntrinsicOrd + AddAssign {
}
}
impl Real for f32 {
/// This exists to allow configuration over # of bits (f32 v.s. f64) as well as
/// constraints.
pub trait Real: ToFloat + decorum_Real + IntrinsicOrd + AddAssign + MulAssign + fmt::LowerExp + fmt::Display {
fn from_f32(f: f32) -> Self {
Self::from_f64(f as _)
}
fn from_f64(f: f64) -> Self {
Self::from_f32(f as _)
}
fn from_u64(u: u64) -> Self {
Self::from_f64(u as _)
}
fn from_u32(u: u32) -> Self {
Self::from_u64(u as _)
}
fn from_primitive<P: ToFloat>(p: P) -> Self {
// TODO: fold with from_<blah>
Self::from_f64(p.to_f64())
}
fn tenth() -> Self {
Self::one() / Self::ten()
}
fn third() -> Self {
Self::one() / Self::three()
}
fn half() -> Self {
Self::from_f32(0.5)
}
fn two() -> Self {
Self::from_f32(2.0)
}
fn three() -> Self {
Self::from_f32(3.0)
}
fn ten() -> Self {
Self::from_f32(10.0)
}
}
impl ToFloat for f32 {
fn to_f32(&self) -> f32 {
*self
}
}
impl Real for f32 {
fn from_f32(f: f32) -> Self {
f
}
fn from_u64(u: u64) -> Self {
Self::from_f32(u as _)
}
}
impl ToFloat for f64 {
fn to_f64(&self) -> f64 {
*self
}
}
impl Real for f64 {
fn to_f64(&self) -> f64 {
*self
fn from_f64(f: f64) -> Self {
f
}
}
impl ToFloat for decorum::R32 {
fn to_f32(&self) -> f32 {
self.into_inner()
}
}
impl Real for decorum::R32 {
fn to_f32(&self) -> f32 {
fn from_f32(f: f32) -> Self {
Self::from_inner(f)
}
fn from_u64(u: u64) -> Self {
Self::from_f32(u as _)
}
}
impl ToFloat for decorum::R64 {
fn to_f64(&self) -> f64 {
self.into_inner()
}
}
impl Real for decorum::R64 {
fn to_f64(&self) -> f64 {
self.into_inner()
fn from_f64(f: f64) -> Self {
Self::from_inner(f)
}
}
impl Real for decorum::N32 {
impl ToFloat for decorum::N32 {
fn to_f32(&self) -> f32 {
self.into_inner()
}
}
impl Real for decorum::N64 {
impl Real for decorum::N32 {
fn from_f32(f: f32) -> Self {
Self::from_inner(f)
}
fn from_u64(u: u64) -> Self {
Self::from_f32(u as _)
}
}
impl ToFloat for decorum::N64 {
fn to_f64(&self) -> f64 {
self.into_inner()
}
}
impl Real for decorum::N64 {
fn from_f64(f: f64) -> Self {
Self::from_inner(f)
}
}
impl ToFloat for u32 {
fn to_f32(&self) -> f32 {
*self as _
}
fn to_f64(&self) -> f64 {
*self as _
}
}
impl ToFloat for u64 {
fn to_f32(&self) -> f32 {
*self as _
}
fn to_f64(&self) -> f64 {
*self as _
}
}

24
src/render.rs
View File

@@ -1,6 +1,7 @@
use crate::geom::{Index, Meters, Vec2, Vec3, Vec3u};
use crate::{Material as _, MaterialExt as _};
use crate::mat;
use crate::real::{Real as _, ToFloat as _};
use crate::sim::{Cell, GenericSim, StaticSim};
use crate::meas::AbstractMeasurement;
use crossterm::{cursor, QueueableCommand as _};
@@ -152,7 +153,7 @@ impl<'a, S: GenericSim> RenderSteps<'a, S> {
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() as f32 + if cell.mat().is_vacuum() {
cell.mat().conductivity().mag().to_f32() + if cell.mat().is_vacuum() {
0.0
} else {
5.0
@@ -193,12 +194,12 @@ impl<'a, S: GenericSim> RenderSteps<'a, S> {
let y_prop = y_px as f32 / self.im.height() as f32;
let y_m = y_prop * (self.sim.height() as f32 * self.sim.feature_size() as f32);
let z_m = self.z as f32 * self.sim.feature_size() as f32;
self.sim.sample(Meters(Vec3::new(x_m as _, y_m as _, z_m as _)))
self.sim.sample(Meters(Vec3::new(x_m, y_m, z_m)))
}
////////////// Ex/Ey/Bz configuration ////////////
fn render_b_z_field(&mut self, scale: f32) {
self.render_scalar_field(1.0e-4 * scale, true, 1, |cell| cell.b().z() as f32);
self.render_scalar_field(1.0e-4 * scale, true, 1, |cell| cell.b().z().to_f32());
}
fn render_e_xy_field(&mut self, scale: f32) {
self.render_vector_field(Rgb([0xff, 0xff, 0xff]), 100.0 * scale, |cell| cell.e().xy().to_f32());
@@ -209,24 +210,24 @@ impl<'a, S: GenericSim> RenderSteps<'a, S> {
}
////////////// Magnitude configuration /////////////
fn render_b(&mut self, scale: f32) {
self.render_scalar_field(1.0e-3 * scale, false, 1, |cell| cell.b().mag() as f32);
self.render_scalar_field(1.0e-3 * scale, false, 1, |cell| cell.b().mag().to_f32());
}
fn render_current(&mut self, scale: f32) {
self.render_scalar_field(1.0e1 * scale, false, 0, |cell| {
cell.e().elem_mul(cell.mat().conductivity()).mag() as f32
cell.e().elem_mul(cell.mat().conductivity()).mag().to_f32()
});
}
////////////// Bx/By/Ez configuration ////////////
fn render_e_z_field(&mut self, scale: f32) {
self.render_scalar_field(1e4 * scale, true, 1, |cell| cell.e().z() as f32);
self.render_scalar_field(1e4 * scale, true, 1, |cell| cell.e().z().to_f32());
}
fn render_b_xy_field(&mut self, scale: f32) {
self.render_vector_field(Rgb([0xff, 0xff, 0xff]), 1.0e-9 * scale, |cell| cell.b().xy().to_f32());
}
fn render_m(&mut self, scale: f32) {
self.render_scalar_field(1.0e5 * scale, false, 1, |cell| cell.mat().m().mag() as f32);
self.render_scalar_field(1.0e5 * scale, false, 1, |cell| cell.mat().m().mag().to_f32());
self.render_vector_field(Rgb([0xff, 0xff, 0xff]), 1.0e5 * scale, |cell| cell.mat().m().xy().to_f32());
}
@@ -558,7 +559,7 @@ impl<S: GenericSim> Renderer<S> for PlotlyRenderer {
yv.push(y);
zv.push(z);
// opacities.push((cell.e().mag() * 0.1).min(1.0) as f64)
let mat = cell.mat().conductivity().mag() + if cell.mat().is_vacuum() {
let mat = cell.mat().conductivity().mag().to_f32() + if cell.mat().is_vacuum() {
0.0
} else {
5.0
@@ -566,9 +567,9 @@ impl<S: GenericSim> Renderer<S> for PlotlyRenderer {
//let g = scale_unsigned_to_u8(mat, 10.0);
//let r = scale_unsigned_to_u8(cell.mat().m().mag(), 100.0);
//let b = scale_unsigned_to_u8(cell.e().mag(), 1e2);
let r = scale_unsigned_to_u8(cell.mat().m().mag() as f32, 100.0);
let g = scale_unsigned_to_u8(cell.e().mag() as f32, 1e2);
let b = scale_unsigned_to_u8(mat as f32, 10.0);
let r = scale_unsigned_to_u8(cell.mat().m().mag().to_f32(), 100.0);
let g = scale_unsigned_to_u8(cell.e().mag().to_f32(), 1e2);
let b = scale_unsigned_to_u8(mat, 10.0);
let alpha = 1.0;
colors.push(Rgba::new(r, g, b, alpha));
}
@@ -666,6 +667,7 @@ impl<S: GenericSim + Clone + Serialize> Renderer<S> for SerializerRenderer {
}
fn render(&self, state: &S, measurements: &[Box<dyn AbstractMeasurement>], _config: RenderConfig) {
let frame = SerializedFrame {
// N.B.: `to_static` isn't 100% necessary here, but it tends to result in a smaller file.
state: state.to_static(),
measurements: measurements.iter().cloned().collect(),
};

174
src/sim.rs
View File

@@ -1,6 +1,7 @@
use crate::{flt::{Flt, Real}, consts};
use crate::geom::{Coord, Cube, Index, InvertedRegion, Meters, Region, Vec3, Vec3u};
use crate::mat::{self, GenericMaterial, Material, MaterialExt as _};
use crate::real::{self, Real as _};
use crate::stim::AbstractStimulus;
use dyn_clone::{self, DynClone};
use log::trace;
@@ -207,7 +208,7 @@ pub trait GenericSim: Send + Sync + DynClone {
}
fn volume(&self) -> Flt {
let s = self.size().to_meters(self.feature_size());
s.x() * s.y() * s.z()
Flt::from_primitive(s.x() * s.y() * s.z())
}
fn timestep(&self) -> Flt;
fn step_no(&self) -> u64;
@@ -309,7 +310,7 @@ impl<'a> dyn GenericSim + 'a {
}
pub fn current<C: Coord>(&self, c: C) -> Vec3 {
self.get(c).current_density() * (self.feature_size() * self.feature_size())
self.get(c).current_density() * Real::from_primitive(self.feature_size() * self.feature_size())
}
pub fn impulse_e<C: Coord>(&mut self, pos: C, amt: Vec3) {
@@ -454,7 +455,7 @@ impl<M: Material + Clone + Default + Send + Sync + 'static> SimState<M> {
Zip::from(ndarray::indices_of(&self.e)).and(&mut self.e).par_apply(
|(z, y, x), e| {
let pos_meters = Index((x as u32, y as u32, z as u32).into()).to_meters(feature_size);
let value = stim.at(t_sec, pos_meters) * timestep;
let value = stim.at(t_sec, pos_meters) * Real::from_primitive(timestep);
*e += value;
});
trace!("apply_stimulus end");
@@ -740,26 +741,26 @@ impl<M> Cell<M> {
self.state.e()
}
pub fn ex(&self) -> Flt {
self.e().x()
Flt::from_primitive(self.e().x())
}
pub fn ey(&self) -> Flt {
self.e().y()
Flt::from_primitive(self.e().y())
}
pub fn ez(&self) -> Flt {
self.e().z()
Flt::from_primitive(self.e().z())
}
pub fn h(&self) -> Vec3 {
self.state.h()
}
pub fn hx(&self) -> Flt {
self.h().x()
Flt::from_primitive(self.h().x())
}
pub fn hy(&self) -> Flt {
self.h().y()
Flt::from_primitive(self.h().y())
}
pub fn hz(&self) -> Flt {
self.h().z()
Flt::from_primitive(self.h().z())
}
}
@@ -768,13 +769,13 @@ impl<M: Material> Cell<M> {
(self.state.h() + self.mat.m()) * consts::real::MU0()
}
pub fn bx(&self) -> Flt {
self.b().x()
Flt::from_primitive(self.b().x())
}
pub fn by(&self) -> Flt {
self.b().y()
Flt::from_primitive(self.b().y())
}
pub fn bz(&self) -> Flt {
self.b().z()
Flt::from_primitive(self.b().z())
}
pub fn mat(&self) -> &M {
@@ -795,17 +796,17 @@ impl<M: Material> Cell<M> {
/// Solve df/dt in:
/// nabla_g = df_dt_coeff * df/dt + f_coeff * f + f_int_coeff * \int f + f_int_int_coeff * \int \int f
#[allow(unused)]
fn solve_step_diff_eq(
nabla_g: Vec3,
df_dt_coeff: Vec3,
f_coeff: Vec3,
f_int_coeff: Vec3,
f_int_int_coeff: Vec3,
delta_t: Real,
f_prev: Vec3,
f_int_prev: Vec3,
f_int_int_prev: Vec3,
) -> Vec3 {
fn solve_step_diff_eq<R: real::Real>(
nabla_g: Vec3<R>,
df_dt_coeff: Vec3<R>,
f_coeff: Vec3<R>,
f_int_coeff: Vec3<R>,
f_int_int_coeff: Vec3<R>,
delta_t: R,
f_prev: Vec3<R>,
f_int_prev: Vec3<R>,
f_int_int_prev: Vec3<R>,
) -> Vec3<R> {
// ```tex
// $f(t)$ from $(T-\Delta t)$ to $(T+\Delta t)$ is a line. As such:
// $f(t') = f_prev + t'*db/dt$ where $t'$ is the time since the last sample of $f$, i.e. $f_prev$
@@ -813,18 +814,16 @@ fn solve_step_diff_eq(
// $\int f|t' = f_int_prev + t'*f_prev + 1/2 t'^2 df/dt$
// $\int \int f|t' = f_int_int_prev + t'*f_int_prev + 1/2 t'^2*f_prev + 1/6 t'^3*df/dt$
// ```
use consts::real::{HALF, THIRD};
// Solve coeffs in:
// nabla_g = df_dt_coeff * df/dt + f_prev_coeff * f_prev
// + f_int_prev_coeff * f_int_prev + f_int_int_prev_coeff * f_int_int_prev
let f_int_int_prev_coeff = f_int_int_coeff;
let f_int_prev_coeff = f_int_coeff + f_int_int_coeff*delta_t;
let f_prev_coeff = f_coeff + (f_int_coeff +
f_int_int_coeff*HALF()
f_int_int_coeff*R::half()
)*delta_t;
let df_dt_coeff = df_dt_coeff + (f_coeff +
(f_int_coeff + f_int_int_coeff*THIRD()*delta_t)*HALF()*delta_t
(f_int_coeff + f_int_int_coeff*R::third()*delta_t)*R::half()*delta_t
)*delta_t;
// Solve coeffs in nabla_g = df_dt_coeff + offset;
@@ -1161,6 +1160,7 @@ mod test {
use super::*;
use crate::geom::{Cube, Region, WorldRegion};
use crate::mat::{Conductor, Static};
use crate::real::Zero as _;
use float_eq::assert_float_eq;
use more_asserts::*;
use rand::rngs::StdRng;
@@ -1170,7 +1170,7 @@ mod test {
let mut cf = ConvFlt(Real::from(0.0));
for (s, e) in signal.iter().zip(expected.iter()) {
let v = cf.step(a.into(), b.into(), c.into(), s.clone().into(), delta_t.into());
assert_float_eq!(v.into_inner(), e, r2nd <= 1e-6);
assert_float_eq!(Flt::from_primitive(v), e, r2nd <= 1e-6);
}
}
@@ -1273,7 +1273,9 @@ mod test {
do_conv_flt_test(a, b, c, delta_t, &signal, &expected);
}
fn assert_vec3_eq(actual: Vec3, expected: Vec3) {
fn assert_vec3_eq<R: real::Real>(actual: Vec3<R>, expected: Vec3<R>) {
let actual = actual.cast::<f64>();
let expected = expected.cast::<f64>();
assert_float_eq!(actual.x(), expected.x(), r2nd <= 1e-6);
assert_float_eq!(actual.y(), expected.y(), r2nd <= 1e-6);
assert_float_eq!(actual.z(), expected.z(), r2nd <= 1e-6);
@@ -1285,11 +1287,11 @@ mod test {
// nabla_g = 0.1 df/dt
let df_dt = solve_step_diff_eq(
Vec3::new(2.0, 3.0, 4.0), // nabla_g
Vec3::unit()*0.1, // df_dt_coeff
Vec3::uniform(0.1), // df_dt_coeff
Vec3::zero(),
Vec3::zero(),
Vec3::zero(),
consts::real::HALF(), // delta_t
Real::half(), // delta_t
Vec3::zero(),
Vec3::zero(),
Vec3::zero(),
@@ -1303,11 +1305,11 @@ mod test {
// nabla_g = 0.1 df/dt - f - 3 f_int + 0.4 f_int_int
let df_dt = solve_step_diff_eq(
Vec3::new(2.0, 3.0, 4.0), // nabla_g
Vec3::unit()*0.1, // df_dt_coeff
Vec3::unit()*-1.0, // f_coeff
Vec3::unit()*-3.0, // f_int_coeff
Vec3::unit()*0.4, // f_int_int_coeff
consts::real::ZERO(), // delta_t
Vec3::uniform(0.1), // df_dt_coeff
Vec3::uniform(-1.0), // f_coeff
Vec3::uniform(-3.0), // f_int_coeff
Vec3::uniform(0.4), // f_int_int_coeff
Real::zero(), // delta_t
Vec3::zero(),
Vec3::zero(),
Vec3::zero(),
@@ -1320,16 +1322,16 @@ mod test {
// nabla_g = 0.1 df/dt + 2 f
// Let f(t) = sin(2 t)
// then f'(t) = 2 cos(2 t)
let f = Vec3::unit()*(5.0).sin();
let df_dt = Vec3::unit()*2.0 * (5.0).cos();
let nabla_g = df_dt*0.1 + f*2.0;
let f = Vec3::uniform((5.0).sin());
let df_dt = Vec3::uniform(2.0 * (5.0).cos());
let nabla_g = df_dt*Real::tenth() + f*Real::two();
let actual_df_dt = solve_step_diff_eq(
nabla_g,
Vec3::unit()*0.1,
Vec3::unit()*2.0,
Vec3::uniform(0.1),
Vec3::uniform(2.0),
Vec3::zero(),
Vec3::zero(),
consts::real::ZERO(), // delta_t
Real::zero(), // delta_t
f, // f_prev
Vec3::zero(),
Vec3::zero(),
@@ -1343,17 +1345,17 @@ mod test {
// Let f(t) = sin(2 t)
// then f'(t) = 2 cos(2 t)
// then \int f(t) = -1/2 cos(2 t)
let f = Vec3::unit()*(5.0).sin();
let df_dt = Vec3::unit()*2.0 * (5.0).cos();
let f_int = Vec3::unit()*-0.5 * (5.0).cos();
let nabla_g = df_dt*0.1 + f*2.0 + f_int*0.4;
let f = Vec3::uniform((5.0).sin());
let df_dt = Vec3::uniform(2.0 * (5.0).cos());
let f_int = Vec3::uniform(-0.5 * (5.0).cos());
let nabla_g = df_dt*Real::tenth() + f*Real::two() + f_int*Real::from_primitive(0.4);
let actual_df_dt = solve_step_diff_eq(
nabla_g,
Vec3::unit()*0.1,
Vec3::unit()*2.0,
Vec3::unit()*0.4,
Vec3::uniform(0.1),
Vec3::uniform(2.0),
Vec3::uniform(0.4),
Vec3::zero(),
consts::real::ZERO(), // delta_t
Real::zero(), // delta_t
f, // f_prev
f_int, // f_int_prev
Vec3::zero(),
@@ -1375,11 +1377,11 @@ mod test {
let nabla_g = df_dt*0.1 + f*2.0 + f_int*0.4 + f_int_int*0.3;
let actual_df_dt = solve_step_diff_eq(
nabla_g,
Vec3::unit()*0.1,
Vec3::unit()*2.0,
Vec3::unit()*0.4,
Vec3::unit()*0.3,
consts::real::ZERO(), // delta_t
Vec3::uniform(0.1),
Vec3::uniform(2.0),
Vec3::uniform(0.4),
Vec3::uniform(0.3),
f32::zero(), // delta_t
f, // f_prev
f_int, // f_int_prev
f_int_int, // f_int_int_prev
@@ -1391,16 +1393,16 @@ mod test {
fn step_diff_eq_1st_order_homogenous() {
// 0 = f'(t) + f(t)
// solution: f(t) = k*e^-t
let f = |t: Flt| -> Vec3 {
let f = |t: f32| -> Vec3<f32> {
Vec3::unit()*(-t).exp()
};
let df_dt = |t: Flt| -> Vec3 {
let df_dt = |t: f32| -> Vec3<f32> {
Vec3::unit()*-(-t).exp()
};
let f_int = |t: Flt| -> Vec3 {
let f_int = |t: f32| -> Vec3<f32> {
Vec3::unit()*-(-t).exp()
};
let f_int_int = |t: Flt| -> Vec3 {
let f_int_int = |t: f32| -> Vec3<f32> {
Vec3::unit()*(-t).exp()
};
let tl = 3.0000;
@@ -1411,7 +1413,7 @@ mod test {
Vec3::unit(), // f_coeff
Vec3::zero(),
Vec3::zero(),
(tm - tl).into(), // delta_t
tm - tl, // delta_t
f(tl), // f_prev
f_int(tl),
f_int_int(tl),
@@ -1423,16 +1425,16 @@ mod test {
fn step_diff_eq_1st_order_nonhomogenous() {
// 1 = f'(t) + f(t)
// let f(t) = 1 + e^-t
let f = |t: Flt| -> Vec3 {
let f = |t: f32| -> Vec3<f32> {
Vec3::unit() + Vec3::unit()*(-t).exp()
};
let df_dt = |t: Flt| -> Vec3 {
let df_dt = |t: f32| -> Vec3<f32> {
Vec3::unit()*-(-t).exp()
};
let f_int = |t: Flt| -> Vec3 {
let f_int = |t: f32| -> Vec3<f32> {
Vec3::unit()*-(-t).exp()
};
let f_int_int = |t: Flt| -> Vec3 {
let f_int_int = |t: f32| -> Vec3<f32> {
Vec3::unit()*(-t).exp()
};
let tl = 3.0000;
@@ -1443,7 +1445,7 @@ mod test {
Vec3::unit(), // f_coeff
Vec3::zero(),
Vec3::zero(),
(tm - tl).into(), // delta_t
tm - tl, // delta_t
f(tl), // f_prev
f_int(tl),
f_int_int(tl),
@@ -1464,16 +1466,16 @@ mod test {
let b = Vec3::new(0.5, -0.3, 0.0);
let c = Vec3::new(-5.0, 0.0, -0.1);
let d = -w.elem_mul(c + w.elem_mul(b + w.elem_mul(a)));
let f = |t: Flt| -> Vec3 {
let f = |t: f64| -> Vec3<f64> {
V.elem_mul((w*t).exp())
};
let df_dt = |t: Flt| -> Vec3 {
let df_dt = |t: f64| -> Vec3<f64> {
V.elem_mul((w*t).exp()).elem_mul(w)
};
let f_int = |t: Flt| -> Vec3 {
let f_int = |t: f64| -> Vec3<f64> {
V.elem_mul((w*t).exp()).elem_div(w)
};
let f_int_int = |t: Flt| -> Vec3 {
let f_int_int = |t: f64| -> Vec3<f64> {
V.elem_mul((w*t).exp()).elem_div(w).elem_div(w) + C.elem_div(d)
};
let tl = 3.00000000;
@@ -1484,7 +1486,7 @@ mod test {
b, // f_coeff
c, // int_f_coeff
d, // int_int_f_coeff
(tm - tl).into(), // delta_t
tm - tl, // delta_t
f(tl), // f_prev
f_int(tl), // f_int_prev
f_int_int(tl), // f_int_int_prev
@@ -1495,7 +1497,7 @@ mod test {
fn energy(s: &dyn GenericSim) -> Flt {
0.5 * s.feature_volume() * s.map_sum_enumerated(|_pos: Index, cell| {
// println!("{}: {}", _pos, cell.e());
Real::from(cell.e().mag_sq()).into_inner()
Flt::from_primitive(cell.e().mag_sq())
})
}
@@ -1582,7 +1584,7 @@ mod test {
#[test]
fn anisotropic_conductor_inactive_x() {
let (energy_0, energy_1) = conductor_test(Conductor::new_anisotropic(
(1e3, 0e3, 0e3).into()
Vec3::new_x(1e3)
));
assert_gt!(energy_1, 0.9*energy_0);
// XXX: I'm not sure why it _gains_ energy :o
@@ -1592,7 +1594,7 @@ mod test {
#[test]
fn anisotropic_conductor_inactive_y() {
let (energy_0, energy_1) = conductor_test(Conductor::new_anisotropic(
(0e3, 1e3, 0e3).into()
Vec3::new_y(1e3)
));
assert_gt!(energy_1, 0.9*energy_0);
assert_lt!(energy_1, 1.5*energy_0);
@@ -1601,7 +1603,7 @@ mod test {
#[test]
fn anisotropic_conductor_active_z() {
let (energy_0, energy_1) = conductor_test(Conductor::new_anisotropic(
(0e3, 0e3, 1e3).into()
Vec3::new_z(1e3)
));
assert_float_eq!(energy_1/energy_0, 0.0, abs <= 1e-6);
}
@@ -1610,8 +1612,8 @@ mod test {
let mut state = StaticSim::new(size, 1e-6);
let timestep = state.timestep();
state.fill_boundary_using(size/4, |boundary_ness| {
let b = boundary_ness.elem_pow(3.0);
let conductivity = b * (0.5 / timestep);
let b = boundary_ness.elem_pow(Real::three());
let conductivity = b * (Real::half() / timestep);
// K scales up to 11; dc_shift scales from 0.05:
// https://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.461.4606&rep=rep1&type=pdf
// let propagation = Vec3::unit() + b * 10.0;
@@ -1629,8 +1631,8 @@ mod test {
let mut state = StaticSim::new(size, 1e-6);
let timestep = state.timestep();
state.fill_boundary_using(size/4, |boundary_ness| {
let b = boundary_ness.elem_pow(3.0);
let conductivity = Vec3::unit() * (consts::real::EPS0() * b.mag() * 0.5 / timestep);
let b = boundary_ness.elem_pow(Real::three());
let conductivity = Vec3::unit() * (consts::real::EPS0() * b.mag() * Real::half() / timestep);
Static {
conductivity, ..Default::default()
}
@@ -1652,7 +1654,7 @@ mod test {
let (e, hz) = (self.f)(pos.to_index(self.feat_size));
let sig_angle = angle*hz;
let sig = sig_angle.sin();
e*gate*sig / (self.t_end * self.sqrt_vol)
e*Real::from_primitive(gate*sig / (self.t_end * self.sqrt_vol))
}
}
@@ -1774,14 +1776,14 @@ mod test {
let mut state = StaticSim::new(size, 1e-6);
let timestep = state.timestep();
state.fill_boundary_using(size/4, |boundary_ness| {
let b = boundary_ness.elem_pow(3.0);
let coord_stretch = b * (0.5 / timestep);
let b = boundary_ness.elem_pow(Real::three());
let coord_stretch = b * (Real::half() / timestep);
// let coord_stretch = b * 0.5;
// Force the stretch to be only along one axis
let coord_stretch = match coord_stretch.to_tuple() {
(x, y, z) if x >= y && x >= z => (x, 0.0, 0.0),
(x, y, z) if y >= x && y >= z => (0.0, y, 0.0),
(x, y, z) if z >= x && z >= y => (0.0, 0.0, z),
(x, y, z) if x >= y && x >= z => (x, Real::zero(), Real::zero()),
(x, y, z) if y >= x && y >= z => (Real::zero(), y, Real::zero()),
(x, y, z) if z >= x && z >= y => (Real::zero(), Real::zero(), z),
_ => unreachable!(),
}.into();
Static::from_pml(coord_stretch)
@@ -1796,8 +1798,8 @@ mod test {
let mut state = StaticSim::new(size, 1e-6);
let timestep = state.timestep();
state.fill_boundary_using(size/4, |boundary_ness| {
let b = boundary_ness.elem_pow(3.0);
let cs = b * (0.5 / timestep);
let b = boundary_ness.elem_pow(Real::three());
let cs = b * (Real::half() / timestep);
// let cs = b * 0.5;
// permute the stretching so that it shouldn't interfere with the wave
let coord_stretch = shuffle(cs);

19
src/stim.rs
View File

@@ -1,5 +1,6 @@
use crate::consts;
use crate::flt::Flt;
use crate::flt::{self, Flt};
use crate::real::*;
use crate::geom::{Meters, Region, Vec3};
pub trait AbstractStimulus: Sync {
@@ -66,7 +67,7 @@ impl<R: Region + Sync, T: TimeVarying1 + Sync> AbstractStimulus for CurlStimulus
let amount = self.stim.at(t_sec);
let from_center_to_point = *pos - *self.center;
let rotational = from_center_to_point.cross(*self.axis);
let impulse = rotational.with_mag(amount);
let impulse = rotational.with_mag(flt::Real::from_primitive(amount));
impulse
} else {
Vec3::zero()
@@ -140,7 +141,7 @@ impl TimeVarying1 for Sinusoid1 {
impl TimeVarying3 for Sinusoid3 {
fn at(&self, t_sec: Flt) -> Vec3 {
self.amp * (t_sec * self.omega).sin()
self.amp * flt::Real::from_primitive(t_sec * self.omega).sin()
}
}
@@ -216,19 +217,19 @@ mod test {
#[test]
fn sinusoid3() {
let s = Sinusoid3::new((10.0, 1.0, -100.0).into(), 1000.0);
let s = Sinusoid3::new(Vec3::new(10.0, 1.0, -100.0), 1000.0);
assert_eq!(s.at(0.0), Vec3::zero());
assert_approx_eq!(s.at(0.00025), (10.0, 1.0, -100.0).into());
assert_approx_eq!(s.at(0.00025), Vec3::new(10.0, 1.0, -100.0));
assert_approx_eq!(s.at(0.00050), Vec3::zero());
assert_approx_eq!(s.at(0.00075), (-10.0, -1.0, 100.0).into());
assert_approx_eq!(s.at(0.00075), Vec3::new(-10.0, -1.0, 100.0));
}
#[test]
fn sinusoid3_from_wavelength() {
let s = Sinusoid3::from_wavelength((10.0, 1.0, -100.0).into(), 0.001);
let s = Sinusoid3::from_wavelength(Vec3::new(10.0, 1.0, -100.0), 0.001);
assert_eq!(s.at(0.0), Vec3::zero());
assert_approx_eq!(s.at(0.00025), (10.0, 1.0, -100.0).into());
assert_approx_eq!(s.at(0.00025), Vec3::new(10.0, 1.0, -100.0));
assert_approx_eq!(s.at(0.00050), Vec3::zero());
assert_approx_eq!(s.at(0.00075), (-10.0, -1.0, 100.0).into());
assert_approx_eq!(s.at(0.00075), Vec3::new(-10.0, -1.0, 100.0));
}
}