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

Stabilize time as f32.

Browse Source
This commit is contained in:
Colin
2021-06-07 18:25:11 -07:00
parent 901d2334f1
commit 8779371f42
12 changed files with 288 additions and 261 deletions
Download Patch File Download Diff File Expand all files Collapse all files

16
examples/minimal_torus.rs
View File

@@ -4,7 +4,7 @@ use coremem::stim::{CurlStimulus, Sinusoid1, TimeVarying1 as _};
fn main() {
coremem::init_logging();
let feat_size = 10e-6; // feature size
let feat_size = 10e-6f32; // feature size
let duration = 6e-9;
let width = 2200e-6;
let depth = 1800e-6;
@@ -16,10 +16,10 @@ fn main() {
let peak_current = 5e11;
let current_duration = 1.0e-11; // half-wavelength of the sine wave
let current_break = 0.5e-11; // time between 'set' pulse and 'clear' pulse
let conductivity = 1.0e9;
let conductivity = 1.0e9f32;
let from_m = |m: Flt| (m/feat_size).round() as u32;
let m_to_um = |m: Flt| (m * 1e6).round() as u32;
let from_m = |m: f32| (m/feat_size).round() as u32;
let m_to_um = |m: f32| (m * 1e6).round() as u32;
let half_width = width * 0.5;
let half_depth = depth * 0.5;
@@ -49,8 +49,8 @@ fn main() {
);
driver.fill_region(&ferro_region, mat::db::minimal_square_ferrite().into());
driver.fill_region(&drive_region, mat::db::conductor(conductivity).into());
driver.fill_region(&sense_region, mat::db::conductor(conductivity).into());
driver.fill_region(&drive_region, mat::db::conductor(conductivity as _).into());
driver.fill_region(&sense_region, mat::db::conductor(conductivity as _).into());
// driver.add_pml_boundary(Meters((boundary_xy, boundary_xy, boundary_z).into()));
driver.add_classical_boundary(Meters::new(boundary_xy, boundary_xy, boundary_z));
@@ -61,10 +61,10 @@ fn main() {
// 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 / (drive_region.cross_section() * conductivity);
let pos_wave = Sinusoid1::from_wavelength(peak_stim, current_duration * 2.0)
let pos_wave = Sinusoid1::from_wavelength(peak_stim as _, current_duration * 2.0)
.half_cycle();
let neg_wave = Sinusoid1::from_wavelength(-peak_stim, current_duration * 2.0)
let neg_wave = Sinusoid1::from_wavelength(-peak_stim as _, current_duration * 2.0)
.half_cycle()
.shifted(current_duration + current_break);

12
examples/toroid25d.rs
View File

@@ -6,11 +6,11 @@ use log::trace;
fn main() {
coremem::init_logging();
for p in 0..20 {
let ferro_depth = 10e-6 * (20-p) as Flt;
let feat_size = 10e-6; // feature size
let ferro_depth = 10e-6 * (20-p) as f32;
let feat_size = 10e-6f32; // feature size
let from_m = |m| (m/feat_size) as u32;
let m_to_um = |px| (px * 1e6) as u32;
let to_m = |px| px as Flt * feat_size;
let to_m = |px| px as f32 * feat_size;
let width = 2500e-6;
let depth = 200e-6;
let conductor_inner_rad = 0e-6;
@@ -91,7 +91,7 @@ fn main() {
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);
let center = Vec2::new(half_width, half_width);
for z_px in 0..depth_px {
for y_px in 0..width_px {
@@ -99,9 +99,9 @@ fn main() {
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 Flt..conductor_outer_rad as Flt).contains(&r) {
if (conductor_inner_rad..conductor_outer_rad).contains(&r) {
driver.state.put_material(loc, mat::db::conductor(conductivity).into());
} else if (ferro_inner_rad as Flt..ferro_outer_rad as Flt).contains(&r) {
} else if (ferro_inner_rad..ferro_outer_rad).contains(&r) {
let half_depth_px = from_m(half_depth);
let ferro_depth_px = from_m(ferro_depth);
if (half_depth_px-ferro_depth_px/2 .. half_depth_px+(ferro_depth_px+1)/2).contains(&z_px) {

22
examples/wrapped_torus.rs
View File

@@ -4,7 +4,7 @@ use coremem::stim::{CurlStimulus, Sinusoid1, TimeVarying1 as _};
fn main() {
coremem::init_logging();
let feat_size = 10e-6; // feature size
let feat_size = 10e-6f32; // feature size
let duration = 0.5e-9;
let width = 4800e-6;
let height = 2200e-6;
@@ -17,11 +17,11 @@ 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 drive_conductivity = 5.0;
let sense_conductivity = 5.0;
let drive_conductivity = 5.0f32;
let sense_conductivity = 5.0f32;
let from_m = |m: Flt| (m/feat_size).round() as u32;
let m_to_um = |m: Flt| (m * 1e6).round() as u32;
let from_m = |m: f32| (m/feat_size).round() as u32;
let m_to_um = |m: f32| (m * 1e6).round() as u32;
let half_width = width * 0.5;
let q1_width = width * 0.25;
let q3_width = width * 0.75;
@@ -42,16 +42,16 @@ fn main() {
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(drive_conductivity).into());
driver.fill_region(&sense1_region, mat::db::conductor(sense_conductivity).into());
driver.fill_region(&drive1_region, mat::db::conductor(drive_conductivity as _).into());
driver.fill_region(&sense1_region, mat::db::conductor(sense_conductivity as _).into());
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(drive_conductivity).into());
driver.fill_region(&sense2_region, mat::db::conductor(sense_conductivity).into());
driver.fill_region(&drive2_region, mat::db::conductor(drive_conductivity as _).into());
driver.fill_region(&sense2_region, mat::db::conductor(sense_conductivity as _).into());
let boundary_xy = q1_width - ferro_major - ferro_minor - buffer;
let boundary_z = half_depth - wire_major - wire_minor - buffer;
@@ -65,10 +65,10 @@ fn main() {
// 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() * drive_conductivity);
let pos_wave = Sinusoid1::from_wavelength(peak_stim, current_duration * 2.0)
let pos_wave = Sinusoid1::from_wavelength(peak_stim as Flt, current_duration * 2.0)
.half_cycle();
let neg_wave = Sinusoid1::from_wavelength(-4.0*peak_stim, current_duration * 2.0)
let neg_wave = Sinusoid1::from_wavelength(-4.0*peak_stim as Flt, current_duration * 2.0)
.half_cycle()
.shifted(current_duration + current_break);

39
src/bin/pml_tuning.rs
View File

@@ -1,18 +1,18 @@
use coremem::*;
use coremem::geom::*;
use coremem::real::Real as _;
use coremem::real::{Real as _, ToFloat 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 {
let e = Real::half() * s.map_sum_over_enumerated(reg, |_pos: Index, cell| {
cell.e().mag_sq()
fn energy<R: Region>(s: &dyn GenericSim, reg: &R) -> f32 {
let e = f64::half() * s.map_sum_over_enumerated(reg, |_pos: Index, cell| {
cell.e().mag_sq().to_f64()
});
Flt::from_primitive(e)
e.cast()
}
fn energy_now_and_then<R: Region>(state: &mut StaticSim, reg: &R, frames: u32) -> (Flt, Flt) {
fn energy_now_and_then<R: Region>(state: &mut StaticSim, reg: &R, frames: u32) -> (f32, f32) {
let energy_0 = energy(state, reg);
for _ in 0..frames {
state.step();
@@ -22,13 +22,14 @@ fn energy_now_and_then<R: Region>(state: &mut StaticSim, reg: &R, frames: u32) -
}
struct PmlStim<F> {
/// Maps index -> (stim vector, stim frequency)
f: F,
t_end: Flt,
feat_size: Flt,
t_end: f32,
feat_size: f32,
}
impl<F: Fn(Index) -> (Vec3, Flt) + Sync> AbstractStimulus for PmlStim<F> {
fn at(&self, t_sec: Flt, pos: Meters) -> Vec3 {
let angle = (t_sec as Flt)/self.t_end*2.0*consts::PI;
impl<F: Fn(Index) -> (Vec3, f32) + Sync> AbstractStimulus for PmlStim<F> {
fn at(&self, t_sec: f32, pos: Meters) -> Vec3 {
let angle = t_sec/self.t_end*f32::two_pi();
let gate = 0.5*(1.0 - angle.cos());
let (e, hz) = (self.f)(pos.to_index(self.feat_size));
let sig_angle = angle*hz;
@@ -39,11 +40,11 @@ impl<F: Fn(Index) -> (Vec3, Flt) + Sync> AbstractStimulus for PmlStim<F> {
/// Apply some stimulus, and then let it decay and measure the ratio of energy left in the system
fn apply_stim_full_interior<F>(state: &mut StaticSim, frames: u32, f: F)
where F: Fn(Index) -> (Vec3, Flt) + Sync // returns (E vector, omega)
where F: Fn(Index) -> (Vec3, f32) + Sync // returns (E vector, omega)
{
let stim = PmlStim {
f,
t_end: (frames as Flt) * state.timestep(),
t_end: (frames as f32) * state.timestep(),
feat_size: state.feature_size(),
};
for _t in 0..frames {
@@ -55,7 +56,7 @@ where F: Fn(Index) -> (Vec3, Flt) + Sync // returns (E vector, omega)
fn apply_stim_over_region<R, F>(state: &mut StaticSim, frames: u32, reg: R, f: F)
where
R: Region,
F: Fn(Index) -> (Vec3, Flt) + Sync,
F: Fn(Index) -> (Vec3, f32) + Sync,
{
let feat = state.feature_size();
apply_stim_full_interior(state, frames, |idx| {
@@ -84,7 +85,7 @@ fn apply_chaotic_stim_over_region<R: Region>(state: &mut StaticSim, frames: u32,
})
}
fn chaotic_pml_test(state: &mut StaticSim, boundary: u32, padding: u32, frames: u32) -> Flt {
fn chaotic_pml_test(state: &mut StaticSim, boundary: u32, padding: u32, frames: u32) -> f32 {
let feat = state.feature_size();
{
let upper_left_idx = Index::unit()*padding;
@@ -102,12 +103,12 @@ fn chaotic_pml_test(state: &mut StaticSim, boundary: u32, padding: u32, frames:
energy_1/energy_0
}
fn state_for_pml(size: Index, boundary: Index, feat_size: Flt, sc_coeff: Flt, cond_coeff: Flt, pow: Flt) -> StaticSim {
fn state_for_pml(size: Index, boundary: Index, feat_size: f32, sc_coeff: f32, cond_coeff: f32, pow: f32) -> StaticSim {
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(Real::from_primitive(pow));
let coord_stretch = b * Real::from_primitive(sc_coeff / timestep);
let b = boundary_ness.elem_pow(pow);
let coord_stretch = b * sc_coeff / timestep;
let conductivity = Vec3::unit() * (b.mag() * cond_coeff / timestep);
unimplemented!()
// Static {
@@ -168,7 +169,7 @@ fn main() {
0.5,
] {
//for cond_coeff in vec![0.0, 1.0, 1e3, 1e6, 1e9] {
for cond_coeff in vec![0.0, 0.5*consts::EPS0] {
for cond_coeff in vec![0.0, 0.5*f32::eps0()] {
for frames in vec![400, 1200] {
for pml_width in vec![1, 2, 4, 8, 16] {
for feat_size in vec![1e-6] {

18
src/driver.rs
View File

@@ -34,7 +34,7 @@ pub struct Driver<M=GenericMaterial> {
}
impl<M: Default> Driver<M> {
pub fn new<C: Coord>(size: C, feature_size: Flt) -> Self {
pub fn new<C: Coord>(size: C, feature_size: f32) -> Self {
Driver {
state: SimState::new(size.to_index(feature_size), feature_size),
renderer: Arc::new(MultiRenderer::new()),
@@ -180,7 +180,7 @@ impl<M: Material + Clone + Default + Send + Sync + 'static> Driver<M> {
}
}
pub fn step_until(&mut self, deadline: Flt) {
pub fn step_until(&mut self, deadline: f32) {
while self.dyn_state().time() < deadline {
self.step();
}
@@ -194,9 +194,9 @@ 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(flt::Real::three());
let conductivity = b * (flt::Real::half() / timestep);
Pml::new(conductivity).into()
let b = boundary_ness.elem_pow(3.0);
let conductivity = b * (0.5 / timestep);
Pml::new(conductivity.cast()).into()
});
}
}
@@ -205,8 +205,8 @@ 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(flt::Real::three());
let cond = b * (flt::Real::half() / timestep);
let b = boundary_ness.elem_pow(3.0);
let cond = b * (0.5 / timestep);
let iso_cond = cond.x() + cond.y() + cond.z();
let iso_conductor = mat::db::conductor(Flt::from_primitive(iso_cond));
@@ -223,8 +223,8 @@ struct StimuliAdapter {
}
impl AbstractStimulus for StimuliAdapter {
fn at(&self, t_sec: Flt, pos: Meters) -> Vec3 {
self.stim.at(t_sec, pos) * flt::Real::from_f64(self.frame_interval as f64)
fn at(&self, t_sec: f32, pos: Meters) -> Vec3 {
self.stim.at(t_sec, pos) * flt::Real::from_primitive(self.frame_interval)
}
}

56
src/geom/region.rs
View File

@@ -1,4 +1,3 @@
use crate::flt::{Flt, Real};
use crate::geom::{Meters, Vec2, Vec3};
use crate::real::Real as _;
@@ -16,16 +15,16 @@ dyn_clone::clone_trait_object!(Region);
#[derive(Copy, Clone, Serialize, Deserialize)]
pub struct CylinderZ {
center: Vec2,
radius: Real,
center: Vec2<f32>,
radius: f32,
}
impl CylinderZ {
// TODO: should be unit-typed
pub fn new(center: Vec2, radius: Flt) -> Self {
pub fn new(center: Vec2<f32>, radius: f32) -> Self {
Self {
center,
radius: radius.into()
radius,
}
}
}
@@ -49,24 +48,24 @@ pub struct Torus {
/// Unit-length vector describing the axis of the torus
normal: Meters,
/// Distance from origin to the "center" of the solid part of the torus
major_rad: Real,
major_rad: f32,
/// Distance from a center-point of the torus to its surface
minor_rad: Real,
minor_rad: f32,
}
impl Torus {
pub fn new(center: Meters, normal: Meters, major_rad: Flt, minor_rad: Flt) -> Self {
pub fn new(center: Meters, normal: Meters, major_rad: f32, minor_rad: f32) -> Self {
Self {
center,
normal,
major_rad: Real::from_inner(major_rad),
minor_rad: Real::from_inner(minor_rad),
major_rad,
minor_rad,
}
}
pub fn new_xy(center: Meters, major_rad: Flt, minor_rad: Flt) -> Self {
pub fn new_xy(center: Meters, major_rad: f32, minor_rad: f32) -> Self {
Self::new(center, Meters(Vec3::new(0.0, 0.0, 1.0)), major_rad, minor_rad)
}
pub fn new_xz(center: Meters, major_rad: Flt, minor_rad: Flt) -> Self {
pub fn new_xz(center: Meters, major_rad: f32, minor_rad: f32) -> Self {
Self::new(center, Meters(Vec3::new(0.0, 1.0, 0.0)), major_rad, minor_rad)
}
pub fn center(&self) -> Meters {
@@ -75,9 +74,8 @@ impl Torus {
pub fn axis(&self) -> Meters {
self.normal
}
pub fn cross_section(&self) -> Flt {
use crate::consts::PI;
(self.minor_rad * self.minor_rad * PI).into()
pub fn cross_section(&self) -> f32 {
self.minor_rad * self.minor_rad * f32::pi()
}
}
@@ -102,21 +100,21 @@ impl Region for Torus {
p_on_plane.with_mag(self.major_rad)
};
let distance_to_circle = (rel_p - q).mag();
distance_to_circle < self.minor_rad.into_inner()
distance_to_circle < self.minor_rad
}
}
#[derive(Copy, Clone, Serialize, Deserialize)]
pub struct Sphere {
center: Meters,
rad: Real,
rad: f32,
}
impl Sphere {
pub fn new(center: Meters, rad: Flt) -> Self {
pub fn new(center: Meters, rad: f32) -> Self {
Self {
center,
rad: Real::from_inner(rad),
rad,
}
}
}
@@ -124,7 +122,7 @@ impl Sphere {
#[typetag::serde]
impl Region for Sphere {
fn contains(&self, p: Meters) -> bool {
(*p - *self.center).mag() < self.rad.into_inner()
(*p - *self.center).mag() < self.rad
}
}
@@ -140,23 +138,23 @@ impl Cube {
pub fn new(lower: Meters, upper: Meters) -> Self {
Self { lower, upper }
}
pub fn x_range(&self) -> Range<Flt> {
Flt::from_primitive(self.lower.x())..Flt::from_primitive(self.upper.x())
pub fn x_range(&self) -> Range<f32> {
self.lower.x()..self.upper.x()
}
pub fn y_range(&self) -> Range<Flt> {
Flt::from_primitive(self.lower.y())..Flt::from_primitive(self.upper.y())
pub fn y_range(&self) -> Range<f32> {
self.lower.y()..self.upper.y()
}
pub fn z_range(&self) -> Range<Flt> {
Flt::from_primitive(self.lower.z())..Flt::from_primitive(self.upper.z())
pub fn z_range(&self) -> Range<f32> {
self.lower.z()..self.upper.z()
}
}
#[typetag::serde]
impl Region for Cube {
fn contains(&self, p: Meters) -> bool {
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()))
self.x_range().contains(&p.x()) &&
self.y_range().contains(&p.y()) &&
self.z_range().contains(&p.z())
}
}

37
src/geom/units.rs
View File

@@ -1,20 +1,19 @@
use crate::flt::{Flt, Real};
use crate::real::{Real as _, ToFloat};
use crate::real::ToFloat;
use serde::{Serialize, Deserialize};
use super::{Vec3, Vec3u};
use std::fmt::{self, Display};
use std::ops::{Add, Deref, Div, Mul, Sub};
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;
fn to_meters(&self, feature_size: f32) -> Meters;
fn to_index(&self, feature_size: f32) -> Index;
fn from_meters(other: Meters, feature_size: f32) -> Self;
fn from_index(other: Index, feature_size: f32) -> Self;
fn from_either(i: Index, m: Meters) -> Self;
}
#[derive(Copy, Clone, Debug, Serialize, Deserialize)]
pub struct Meters(pub Vec3);
pub struct Meters(pub Vec3<f32>);
impl Meters {
pub fn new<R: ToFloat>(x: R, y: R, z: R) -> Self {
@@ -23,16 +22,16 @@ impl Meters {
}
impl Coord for Meters {
fn to_meters(&self, _feature_size: Flt) -> Meters {
fn to_meters(&self, _feature_size: f32) -> Meters {
*self
}
fn to_index(&self, feature_size: Flt) -> Index {
Index((self.0 / Real::from_f64(feature_size)).round().into())
fn to_index(&self, feature_size: f32) -> Index {
Index((self.0 / feature_size).round().into())
}
fn from_meters(other: Meters, _feature_size: Flt) -> Self {
fn from_meters(other: Meters, _feature_size: f32) -> Self {
other
}
fn from_index(other: Index, feature_size: Flt) -> Self {
fn from_index(other: Index, feature_size: f32) -> Self {
other.to_meters(feature_size)
}
fn from_either(_i: Index, m: Meters) -> Self {
@@ -41,8 +40,8 @@ impl Coord for Meters {
}
impl Deref for Meters {
type Target = Vec3;
fn deref(&self) -> &Vec3 {
type Target = Vec3<f32>;
fn deref(&self) -> &Self::Target {
&self.0
}
}
@@ -80,16 +79,16 @@ impl Index {
}
impl Coord for Index {
fn to_meters(&self, feature_size: Flt) -> Meters {
Meters(Vec3::from(self.0) * Real::from_f64(feature_size))
fn to_meters(&self, feature_size: f32) -> Meters {
Meters(Vec3::from(self.0) * feature_size)
}
fn to_index(&self, _feature_size: Flt) -> Index {
fn to_index(&self, _feature_size: f32) -> Index {
*self
}
fn from_meters(other: Meters, feature_size: Flt) -> Self {
fn from_meters(other: Meters, feature_size: f32) -> Self {
other.to_index(feature_size)
}
fn from_index(other: Index, _feature_size: Flt) -> Self {
fn from_index(other: Index, _feature_size: f32) -> Self {
other
}
fn from_either(i: Index, _m: Meters) -> Self {

6
src/geom/vecu.rs
View File

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

75
src/meas.rs
View File

@@ -1,6 +1,5 @@
use crate::flt::{Flt, Real};
use crate::geom::{Meters, Region, Torus, Vec3, WorldRegion};
use crate::real::Real as _;
use crate::real::{Real as _, ToFloat as _};
use crate::sim::GenericSim;
use common_macros::b_tree_map;
use dyn_clone::{self, DynClone};
@@ -84,19 +83,19 @@ impl Current {
region: Box::new(r)
}
}
fn data(&self, state: &dyn GenericSim) -> (Real, Vec3) {
fn data(&self, state: &dyn GenericSim) -> (f32, Vec3<f32>) {
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)
FieldSample(1, current.mag().cast(), current.cast())
});
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)
let mean_current_mag = current_mag.to_f32() / (f32::from_primitive(volume));
let mean_current_vec = current_vec.cast::<f32>() / (f32::from_primitive(volume));
(mean_current_mag, mean_current_vec.cast())
}
}
#[derive(Default)]
struct FieldSample(u32, Real, Vec3);
struct FieldSample(u32, f64, Vec3<f64>);
impl std::iter::Sum for FieldSample {
fn sum<I>(iter: I) -> Self
@@ -172,16 +171,16 @@ impl CurrentLoop {
region: r,
}
}
fn data(&self, state: &dyn GenericSim) -> Real {
fn data(&self, state: &dyn GenericSim) -> f32 {
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();
let tangent = normal.cross(to_coord).norm();
let current = state.current(coord);
let directed_current = current.dot(tangent);
FieldSample(1, directed_current, current)
let directed_current = current.dot(tangent.cast());
FieldSample(1, directed_current.cast(), current.cast())
});
let mean_directed_current = directed_current / Real::from_u32(volume);
let mean_directed_current = directed_current.cast::<f32>() / f32::from_primitive(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 +216,7 @@ impl MagneticLoop {
region: r,
}
}
fn data(&self, state: &dyn GenericSim) -> (Real, Real, Real) {
fn data(&self, state: &dyn GenericSim) -> (f32, f32, f32) {
let FieldSamples([
FieldSample(volume, directed_m, _m_vec),
FieldSample(_, directed_b, _b_vec),
@@ -229,23 +228,23 @@ impl MagneticLoop {
let tangent = normal.cross(to_coord).norm();
let m = cell.m();
let directed_m = m.dot(tangent);
let directed_m = m.dot(tangent.cast());
let b = cell.b();
let directed_b = b.dot(tangent);
let directed_b = b.dot(tangent.cast());
let h = cell.h();
let directed_h = h.dot(tangent);
let directed_h = h.dot(tangent.cast());
FieldSamples([
FieldSample(1, directed_m, m),
FieldSample(1, directed_b, b),
FieldSample(1, directed_h, h),
FieldSample(1, directed_m.cast(), m.cast()),
FieldSample(1, directed_b.cast(), b.cast()),
FieldSample(1, directed_h.cast(), h.cast()),
])
});
// let cross_section = self.region.cross_section() / (state.feature_size() * state.feature_size());
let mean_directed_m = directed_m / Real::from_u32(volume);
let mean_directed_m = directed_m.cast::<f32>() / f32::from_primitive(volume);
// let cross_sectional_m = mean_directed_m * cross_section;
let mean_directed_b = directed_b / Real::from_u32(volume);
let mean_directed_b = directed_b.cast::<f32>() / f32::from_primitive(volume);
// let cross_sectional_b = mean_directed_b * cross_section;
let mean_directed_h = directed_h / Real::from_u32(volume);
let mean_directed_h = directed_h.cast::<f32>() / f32::from_primitive(volume);
// let cross_sectional_h = mean_directed_h * cross_section;
// format!(
// "M({}): {:.2e}; B({}): {:.2e}; H({}): {:.2e}",
@@ -292,13 +291,13 @@ impl MagneticFlux {
region: Box::new(r)
}
}
fn data(&self, state: &dyn GenericSim) -> Vec3 {
fn data(&self, state: &dyn GenericSim) -> Vec3<f32> {
let FieldSample(volume, _directed_mag, mag_vec) = state.map_sum_over(&*self.region, |cell| {
let b = cell.b();
let mag = b.mag();
FieldSample(1, mag, b)
FieldSample(1, mag.cast(), b.cast())
});
let mean_mag = mag_vec / Real::from_u32(volume);
let mean_mag = mag_vec.cast() / f32::from_primitive(volume);
mean_mag
}
}
@@ -331,13 +330,13 @@ impl Magnetization {
region: Box::new(r)
}
}
fn data(&self, state: &dyn GenericSim) -> Vec3 {
fn data(&self, state: &dyn GenericSim) -> Vec3<f32> {
let FieldSample(volume, _directed_mag, mag_vec) = state.map_sum_over(&*self.region, |cell| {
let m = cell.m();
let mag = m.mag();
FieldSample(1, mag, m)
FieldSample(1, mag.cast(), m.cast())
});
let mean_mag = mag_vec / Real::from_u32(volume);
let mean_mag = mag_vec.cast() / f32::from_primitive(volume);
mean_mag
}
}
@@ -357,7 +356,7 @@ impl AbstractMeasurement for Magnetization {
}
fn loc(v: Meters) -> String {
format!("{:.0} um", *v * Real::from_u64(1_000_000))
format!("{:.0} um", *v * f32::from_primitive(1_000_000))
}
/// M
@@ -447,7 +446,7 @@ impl Energy {
region: Box::new(region),
}
}
fn data(&self, state: &dyn GenericSim) -> Real {
fn data(&self, state: &dyn GenericSim) -> f32 {
// 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,11 +455,11 @@ impl Energy {
// U(E) = 1/2 \int E . D dV
#[allow(non_snake_case)]
let dV = state.feature_volume();
let e = Real::from_primitive(0.5 * dV) * state.map_sum_over(&*self.region, |cell| {
let e = f64::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()
cell.h().dot(cell.b()).to_f64() + cell.e().mag_sq().to_f64()
});
e
e.cast()
}
}
@@ -494,15 +493,15 @@ impl Power {
region: Box::new(region),
}
}
fn data(&self, state: &dyn GenericSim) -> Flt {
fn data(&self, state: &dyn GenericSim) -> f32 {
// 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 = Real::from_primitive(state.feature_volume());
let power = dV * state.map_sum_over(&*self.region, |cell| {
cell.current_density().dot(cell.e())
let dV = state.feature_volume();
let power = f64::from_primitive(dV) * state.map_sum_over(&*self.region, |cell| {
cell.current_density().dot(cell.e()).to_f64()
});
Flt::from_primitive(power)
power.cast()
}
}

58
src/real.rs
View File

@@ -16,6 +16,11 @@ pub trait ToFloat {
/// 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 {
// TODO: fold with from_<blah>
fn from_primitive<P: ToFloat>(p: P) -> Self {
Self::from_f64(p.to_f64())
}
fn from_f32(f: f32) -> Self {
Self::from_f64(f as _)
}
@@ -23,15 +28,14 @@ pub trait Real: ToFloat + decorum_Real + IntrinsicOrd + AddAssign + MulAssign +
Self::from_f32(f as _)
}
fn from_u64(u: u64) -> Self {
Self::from_f64(u as _)
Self::from_primitive(u)
}
fn from_u32(u: u32) -> Self {
Self::from_u64(u as _)
Self::from_primitive(u)
}
fn from_primitive<P: ToFloat>(p: P) -> Self {
// TODO: fold with from_<blah>
Self::from_f64(p.to_f64())
fn cast<R: Real>(&self) -> R {
R::from_primitive(*self)
}
fn tenth() -> Self {
@@ -52,6 +56,23 @@ pub trait Real: ToFloat + decorum_Real + IntrinsicOrd + AddAssign + MulAssign +
fn ten() -> Self {
Self::from_f32(10.0)
}
fn pi() -> Self {
Self::from_f64(std::f64::consts::PI)
}
fn two_pi() -> Self {
Self::from_f64(std::f64::consts::PI * 2.0)
}
fn c() -> Self {
Self::from_primitive(299792458)
}
fn eps0() -> Self {
Self::from_primitive(8.854187812813e-12) // F⋅m−1
}
fn mu0() -> Self {
Self::from_primitive(1.2566370621219e-6) // H/m
}
}
impl ToFloat for f32 {
@@ -61,12 +82,12 @@ impl ToFloat for f32 {
}
impl Real for f32 {
fn from_primitive<P: ToFloat>(p: P) -> Self {
Self::from_f32(p.to_f32())
}
fn from_f32(f: f32) -> Self {
f
}
fn from_u64(u: u64) -> Self {
Self::from_f32(u as _)
}
}
impl ToFloat for f64 {
@@ -88,12 +109,12 @@ impl ToFloat for decorum::R32 {
}
impl Real for decorum::R32 {
fn from_primitive<P: ToFloat>(p: P) -> Self {
Self::from_f32(p.to_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 {
@@ -115,12 +136,12 @@ impl ToFloat for decorum::N32 {
}
impl Real for decorum::N32 {
fn from_primitive<P: ToFloat>(p: P) -> Self {
Self::from_f32(p.to_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::N64 {
@@ -135,6 +156,15 @@ impl Real for decorum::N64 {
}
}
impl ToFloat for i32 {
fn to_f32(&self) -> f32 {
*self as _
}
fn to_f64(&self) -> f64 {
*self as _
}
}
impl ToFloat for u32 {
fn to_f32(&self) -> f32 {
*self as _

143
src/sim.rs
View File

@@ -1,7 +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::real::{self, Real as _, ToFloat as _};
use crate::stim::AbstractStimulus;
use dyn_clone::{self, DynClone};
use log::trace;
@@ -201,16 +201,16 @@ pub trait GenericSim: Send + Sync + DynClone {
fn impulse_h_meters(&mut self, pos: Meters, amount: Vec3);
fn impulse_b_meters(&mut self, pos: Meters, amount: Vec3);
fn size(&self) -> Index;
fn feature_size(&self) -> Flt;
fn feature_volume(&self) -> Flt {
fn feature_size(&self) -> f32;
fn feature_volume(&self) -> f32 {
let f = self.feature_size();
f*f*f
}
fn volume(&self) -> Flt {
fn volume(&self) -> f32 {
let s = self.size().to_meters(self.feature_size());
Flt::from_primitive(s.x() * s.y() * s.z())
s.x() * s.y() * s.z()
}
fn timestep(&self) -> Flt;
fn timestep(&self) -> f32;
fn step(&mut self);
fn step_no(&self) -> u64;
@@ -223,8 +223,8 @@ pub trait GenericSim: Send + Sync + DynClone {
fn depth(&self) -> u32 {
self.size().z()
}
fn time(&self) -> Flt {
self.timestep() * self.step_no() as Flt
fn time(&self) -> f32 {
self.timestep() * self.step_no() as f32
}
fn apply_stimulus(&mut self, stim: &dyn AbstractStimulus);
@@ -389,16 +389,16 @@ impl<M> SimState<M> {
pub fn depth(&self) -> u32 {
self.size().z()
}
pub fn feature_size(&self) -> Flt {
self.feature_size.into()
pub fn feature_size(&self) -> f32 {
self.feature_size.to_f32()
}
pub fn timestep(&self) -> Flt {
self.timestep.into()
pub fn timestep(&self) -> f32 {
self.timestep.to_f32()
}
}
impl<M: Default> SimState<M> {
pub fn new(size: Index, feature_size: Flt) -> Self {
pub fn new(size: Index, feature_size: f32) -> Self {
// XXX this has to be multiplied by a Courant Number in order to ensure stability.
// For 3d, we need 3 full timesteps in order to communicate information across the corner
// of a grid. The distance traveled during that time is \sqrt{3}. Hence each timestep needs
@@ -406,8 +406,8 @@ impl<M: Default> SimState<M> {
// In 2d, this would be \sqrt{2}/2.
// It's an upper limit though; it should be safe to go lower.
let courant = 0.577;
let feature_size = Real::from(feature_size);
let timestep = feature_size / consts::real::C() * courant;
let feature_size = Real::from_primitive(feature_size);
let timestep = feature_size / Real::c() * Real::from_primitive(courant);
Self {
cells: Array3::default((size.z() as _, size.y() as _, size.x() as _)),
e: Array3::default((size.z() as _, size.y() as _, size.x() as _)),
@@ -456,12 +456,12 @@ impl<M: Material + Clone + Send + Sync + 'static> SimState<M> {
let feature_size = self.feature_size();
let t_sec = self.time();
let timestep = self.timestep();
let timestep = self.timestep;
trace!("apply_stimulus begin");
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) * Real::from_primitive(timestep);
let value = stim.at(t_sec, pos_meters) * timestep;
*e += value;
});
trace!("apply_stimulus end");
@@ -507,7 +507,7 @@ impl<M: Material> SimState<M> {
pub fn fill_boundary_using<C, F>(&mut self, thickness: C, f: F)
where
C: Coord,
F: Fn(Vec3) -> M,
F: Fn(Vec3<f32>) -> M,
{
// TODO: maybe this function belongs on the Driver?
let feat = self.feature_size();
@@ -517,23 +517,23 @@ impl<M: Material> SimState<M> {
let region = InvertedRegion::new(Cube::new(upper_left.to_meters(feat), lower_right.to_meters(feat)));
self.fill_region_using(&region, |loc: Index| {
let depth_x = if loc.x() < upper_left.x() {
(upper_left.x() - loc.x()) as Flt / upper_left.x() as Flt
(upper_left.x() - loc.x()) as f32 / upper_left.x() as f32
} else if loc.x() > lower_right.x() {
(loc.x() - lower_right.x()) as Flt / upper_left.x() as Flt
(loc.x() - lower_right.x()) as f32 / upper_left.x() as f32
} else {
0.0
};
let depth_y = if loc.y() < upper_left.y() {
(upper_left.y() - loc.y()) as Flt / upper_left.y() as Flt
(upper_left.y() - loc.y()) as f32 / upper_left.y() as f32
} else if loc.y() > lower_right.y() {
(loc.y() - lower_right.y()) as Flt / upper_left.y() as Flt
(loc.y() - lower_right.y()) as f32 / upper_left.y() as f32
} else {
0.0
};
let depth_z = if loc.z() < upper_left.z() {
(upper_left.z() - loc.z()) as Flt / upper_left.z() as Flt
(upper_left.z() - loc.z()) as f32 / upper_left.z() as f32
} else if loc.z() > lower_right.z() {
(loc.z() - lower_right.z()) as Flt / upper_left.z() as Flt
(loc.z() - lower_right.z()) as f32 / upper_left.z() as f32
} else {
0.0
};
@@ -584,11 +584,11 @@ impl<M: Material + Clone + Send + Sync + 'static> GenericSim for SimState<M> {
self.cells.shape()[0] as _,
))
}
fn feature_size(&self) -> Flt {
self.feature_size.into()
fn feature_size(&self) -> f32 {
self.feature_size.to_f32()
}
fn timestep(&self) -> Flt {
self.timestep()
fn timestep(&self) -> f32 {
self.timestep.to_f32()
}
fn step(&mut self) {
SimState::step(self)
@@ -601,34 +601,34 @@ impl<M: Material + Clone + Send + Sync + 'static> GenericSim for SimState<M> {
#[allow(unused)]
impl<M> SimState<M> {
fn get_mat<C: Coord>(&self, c: C) -> &M {
let at = c.to_index(self.feature_size.into());
let at = c.to_index(self.feature_size.cast());
&self.cells[at.row_major_idx()]
}
fn get_mat_mut<C: Coord>(&mut self, c: C) -> &mut M {
let at = c.to_index(self.feature_size.into());
let at = c.to_index(self.feature_size.cast());
&mut self.cells[at.row_major_idx()]
}
fn get_e<C: Coord>(&self, c: C) -> &Vec3 {
let at = c.to_index(self.feature_size.into());
let at = c.to_index(self.feature_size.cast());
&self.e[at.row_major_idx()]
}
fn get_e_mut<C: Coord>(&mut self, c: C) -> &mut Vec3 {
let at = c.to_index(self.feature_size.into());
let at = c.to_index(self.feature_size.cast());
&mut self.e[at.row_major_idx()]
}
fn get_h<C: Coord>(&self, c: C) -> &Vec3 {
let at = c.to_index(self.feature_size.into());
let at = c.to_index(self.feature_size.cast());
&self.h[at.row_major_idx()]
}
fn get_h_mut<C: Coord>(&mut self, c: C) -> &mut Vec3 {
let at = c.to_index(self.feature_size.into());
let at = c.to_index(self.feature_size.cast());
&mut self.h[at.row_major_idx()]
}
fn get_hcell<'a, C: Coord>(&'a mut self, c: C) -> HCell<'a, M> {
let at = c.to_index(self.feature_size.into());
let at = c.to_index(self.feature_size.cast());
let idx = at.row_major_idx();
HCell {
e: &self.e[idx],
@@ -1487,14 +1487,14 @@ mod test {
assert_vec3_eq(actual_df_dt, df_dt(tm));
}
fn energy(s: &dyn GenericSim) -> Flt {
fn energy(s: &dyn GenericSim) -> f32 {
0.5 * s.feature_volume() * s.map_sum_enumerated(|_pos: Index, cell| {
// println!("{}: {}", _pos, cell.e());
Flt::from_primitive(cell.e().mag_sq())
})
cell.e().mag_sq().to_f64()
}).to_f32()
}
fn energy_now_and_then<S: GenericSim>(state: &mut S, frames: u32) -> (Flt, Flt) {
fn energy_now_and_then<S: GenericSim>(state: &mut S, frames: u32) -> (f32, f32) {
let energy_0 = energy(state);
for _f in 0..frames {
// println!("e({}) = {}", _f, energy(state));
@@ -1521,7 +1521,7 @@ mod test {
fn energy_conservation_over_time() {
let mut state = StaticSim::new(Index((2001, 1, 1).into()), 1e-6);
for t in 0..100 {
let angle = (t as Flt)/100.0*2.0*consts::PI;
let angle = (t as Flt)/100.0*Flt::two_pi();
let sig = angle.sin();
let gate = 0.5*(1.0 - angle.cos());
//println!("{}", g.exp());
@@ -1539,7 +1539,7 @@ mod test {
fn sane_boundary_conditions() {
let mut state = StaticSim::new(Index((21, 21, 21).into()), 1e-6);
for t in 0..10 {
let angle = (t as Flt)/10.0*2.0*consts::PI;
let angle = (t as Flt)/10.0*Flt::two_pi();
let sig = angle.sin();
let gate = 0.5*(1.0 - angle.cos());
let dyn_state = &mut state as &mut dyn GenericSim;
@@ -1554,11 +1554,11 @@ mod test {
/// Fill the world with the provided material and a stimulus.
/// Measure energy at the start, and then again after advancing many steps.
/// Return these two measurements (energy(t=0), energy(t=~=1000))
fn conductor_test<M: Into<mat::Static>>(mat: M) -> (Flt, Flt) {
fn conductor_test<M: Into<mat::Static>>(mat: M) -> (f32, f32) {
let mut state = StaticSim::new(Index((201, 1, 1).into()), 1e-6);
state.fill_region(&WorldRegion, mat.into());
for t in 0..100 {
let angle = (t as Flt)/100.0*2.0*consts::PI;
let angle = (t as Flt)/100.0*Flt::two_pi();
let sig = angle.sin();
let gate = 0.5*(1.0 - angle.cos());
let dyn_state = &mut state as &mut dyn GenericSim;
@@ -1605,8 +1605,8 @@ mod test {
let mut state = SimState::new(size, 1e-6);
let timestep = state.timestep();
state.fill_boundary_using(size/4, |boundary_ness| {
let b = boundary_ness.elem_pow(Real::three());
let conductivity = b * (Real::half() / timestep);
let b = boundary_ness.elem_pow(3.0);
let conductivity = b * (0.5 / 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;
@@ -1614,7 +1614,7 @@ mod test {
//
// let propagation = Vec3::unit();
// let dc_shift = Vec3::unit() * 1e-6;
Pml::new(conductivity)
Pml::new(conductivity.cast())
});
state
}
@@ -1624,8 +1624,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(Real::three());
let conductivity = Vec3::unit() * (consts::real::EPS0() * b.mag() * Real::half() / timestep);
let b = boundary_ness.elem_pow(3.0);
let conductivity = Vec3::uniform(f32::eps0() * b.mag() * 0.5 / timestep);
Static {
conductivity, ..Default::default()
}
@@ -1635,14 +1635,15 @@ mod test {
/// Stimulus which applies a (possibly different) e*sin(wt)*gate(t) at every point
struct PmlStim<F> {
/// Maps location -> (field direction, frequency)
f: F,
t_end: Flt,
feat_size: Flt,
sqrt_vol: Flt,
t_end: f32,
feat_size: f32,
sqrt_vol: f32,
}
impl<F: Fn(Index) -> (Vec3, Flt) + Sync> AbstractStimulus for PmlStim<F> {
fn at(&self, t_sec: Flt, pos: Meters) -> Vec3 {
let angle = (t_sec as Flt)/self.t_end*2.0*consts::PI;
impl<F: Fn(Index) -> (Vec3, f32) + Sync> AbstractStimulus for PmlStim<F> {
fn at(&self, t_sec: f32, pos: Meters) -> Vec3 {
let angle = t_sec/self.t_end * f32::two_pi();
let gate = 0.5*(1.0 - angle.cos());
let (e, hz) = (self.f)(pos.to_index(self.feat_size));
let sig_angle = angle*hz;
@@ -1652,8 +1653,8 @@ mod test {
}
/// Returns the energy ratio (Eend / Estart)
fn pml_test_full_interior<F, S: GenericSim>(state: &mut S, f: F) -> Flt
where F: Fn(Index) -> (Vec3, Flt) + Sync // returns (E vector, cycles (like Hz))
fn pml_test_full_interior<F, S: GenericSim>(state: &mut S, f: F) -> f32
where F: Fn(Index) -> (Vec3, f32) + Sync // returns (E vector, cycles (like Hz))
{
let stim = PmlStim {
f,
@@ -1670,13 +1671,13 @@ mod test {
// sanity check the starting energy
assert_gt!(energy_0, 1e-11);
assert_lt!(energy_0, 1.0);
(Real::from(energy_1) / Real::from(energy_0)).into_inner()
energy_1 / energy_0
}
fn pml_test_over_region<R, F, S>(state: &mut S, reg: R, f: F) -> Flt
fn pml_test_over_region<R, F, S>(state: &mut S, reg: R, f: F) -> f32
where
R: Region,
F: Fn(Index) -> (Vec3, Flt) + Sync,
F: Fn(Index) -> (Vec3, f32) + Sync,
S: GenericSim
{
let feat = state.feature_size();
@@ -1688,13 +1689,13 @@ mod test {
}
})
}
fn pml_test_uniform_over_region<R: Region, S: GenericSim>(state: &mut S, reg: R, e: Vec3, hz: Flt) -> Flt {
fn pml_test_uniform_over_region<R: Region, S: GenericSim>(state: &mut S, reg: R, e: Vec3, hz: f32) -> f32 {
pml_test_over_region(state, reg, |_idx| {
(e, hz)
})
}
fn pml_test_at<S: GenericSim>(state: &mut S, e: Vec3, center: Index) -> Flt {
fn pml_test_at<S: GenericSim>(state: &mut S, e: Vec3, center: Index) -> f32 {
pml_test_full_interior(state, |idx| {
if idx == center {
(e, 1.0)
@@ -1770,17 +1771,17 @@ mod test {
let mut state = SimState::new(size, 1e-6);
let timestep = state.timestep();
state.fill_boundary_using(size/4, |boundary_ness| {
let b = boundary_ness.elem_pow(Real::three());
let coord_stretch = b * (Real::half() / timestep);
let b = boundary_ness.elem_pow(3.0);
let coord_stretch = b * (0.5 / 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, 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),
let coord_stretch: Vec3<_> = 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),
_ => unreachable!(),
}.into();
Pml::new(coord_stretch)
Pml::new(coord_stretch.cast())
});
state
}
@@ -1792,11 +1793,11 @@ mod test {
let mut state = SimState::new(size, 1e-6);
let timestep = state.timestep();
state.fill_boundary_using(size/4, |boundary_ness| {
let b = boundary_ness.elem_pow(Real::three());
let cs = b * (Real::half() / timestep);
let b = boundary_ness.elem_pow(3.0);
let cs = b * (0.5 / timestep);
// let cs = b * 0.5;
// permute the stretching so that it shouldn't interfere with the wave
let coord_stretch = shuffle(cs);
let coord_stretch = shuffle(cs.cast());
Pml::new(coord_stretch)
});
let center = Index(state.size().0 / 2);

67
src/stim.rs
View File

@@ -1,27 +1,26 @@
use crate::consts;
use crate::flt::{self, Flt};
use crate::real::*;
use crate::geom::{Meters, Region, Vec3};
pub trait AbstractStimulus: Sync {
/// Return the E field which should be added PER-SECOND to the provided position/time.
fn at(&self, t_sec: Flt, pos: Meters) -> Vec3;
fn at(&self, t_sec: f32, pos: Meters) -> Vec3;
}
// impl<T: AbstractStimulus> AbstractStimulus for &T {
// fn at(&self, t_sec: Flt, pos: Meters) -> Vec3 {
// fn at(&self, t_sec: f32, pos: Meters) -> Vec3 {
// (*self).at(t_sec, pos)
// }
// }
impl<T: AbstractStimulus> AbstractStimulus for Vec<T> {
fn at(&self, t_sec: Flt, pos: Meters) -> Vec3 {
fn at(&self, t_sec: f32, pos: Meters) -> Vec3 {
self.iter().map(|s| s.at(t_sec, pos)).sum()
}
}
impl AbstractStimulus for Box<dyn AbstractStimulus> {
fn at(&self, t_sec: Flt, pos: Meters) -> Vec3 {
fn at(&self, t_sec: f32, pos: Meters) -> Vec3 {
(**self).at(t_sec, pos)
}
}
@@ -42,7 +41,7 @@ impl<R, T> Stimulus<R, T> {
}
impl<R: Region + Sync, T: TimeVarying3 + Sync> AbstractStimulus for Stimulus<R, T> {
fn at(&self, t_sec: Flt, pos: Meters) -> Vec3 {
fn at(&self, t_sec: f32, pos: Meters) -> Vec3 {
if self.region.contains(pos) {
self.stim.at(t_sec)
} else {
@@ -68,12 +67,12 @@ impl<R, T> CurlStimulus<R, T> {
}
impl<R: Region + Sync, T: TimeVarying1 + Sync> AbstractStimulus for CurlStimulus<R, T> {
fn at(&self, t_sec: Flt, pos: Meters) -> Vec3 {
fn at(&self, t_sec: f32, pos: Meters) -> Vec3 {
if self.region.contains(pos) {
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(flt::Real::from_primitive(amount));
let rotational: Vec3 = from_center_to_point.cross(*self.axis).cast();
let impulse = rotational.with_mag(amount.cast());
impulse
} else {
Vec3::zero()
@@ -84,22 +83,22 @@ impl<R: Region + Sync, T: TimeVarying1 + Sync> AbstractStimulus for CurlStimulus
pub trait TimeVarying1: Sized {
/// Retrieve the E impulse to apply PER-SECOND at the provided time (in seconds).
fn at(&self, t_sec: Flt) -> Flt;
fn shifted(self, new_start: Flt) -> Shifted<Self> {
fn at(&self, t_sec: f32) -> Flt;
fn shifted(self, new_start: f32) -> Shifted<Self> {
Shifted::new(self, new_start)
}
fn gated(self, from: Flt, to: Flt) -> Gated<Self> {
fn gated(self, from: f32, to: f32) -> Gated<Self> {
Gated::new(self, from, to)
}
}
pub trait TimeVarying3: Sized {
/// Retrieve the E impulse to apply PER-SECOND at the provided time (in seconds).
fn at(&self, t_sec: Flt) -> Vec3;
fn shifted(self, new_start: Flt) -> Shifted<Self> {
fn at(&self, t_sec: f32) -> Vec3;
fn shifted(self, new_start: f32) -> Shifted<Self> {
Shifted::new(self, new_start)
}
fn gated(self, from: Flt, to: Flt) -> Gated<Self> {
fn gated(self, from: f32, to: f32) -> Gated<Self> {
Gated::new(self, from, to)
}
}
@@ -107,26 +106,26 @@ pub trait TimeVarying3: Sized {
#[derive(Clone)]
pub struct Sinusoid<A> {
amp: A,
omega: Flt,
omega: f32,
}
pub type Sinusoid1 = Sinusoid<Flt>;
pub type Sinusoid3 = Sinusoid<Vec3>;
impl<A> Sinusoid<A> {
pub fn new(amp: A, freq: Flt) -> Self {
pub fn new(amp: A, freq: f32) -> Self {
Self {
amp,
omega: freq * consts::TWO_PI,
omega: freq * f32::two_pi(),
}
}
pub fn from_wavelength(amp: A, lambda: Flt) -> Self {
pub fn from_wavelength(amp: A, lambda: f32) -> Self {
Self::new(amp, 1.0/lambda)
}
pub fn freq(&self) -> Flt {
self.omega / consts::TWO_PI
pub fn freq(&self) -> f32 {
self.omega / f32::two_pi()
}
pub fn wavelength(&self) -> Flt {
pub fn wavelength(&self) -> f32 {
1.0 / self.freq()
}
pub fn one_cycle(self) -> Gated<Self> {
@@ -140,13 +139,13 @@ impl<A> Sinusoid<A> {
}
impl TimeVarying1 for Sinusoid1 {
fn at(&self, t_sec: Flt) -> Flt {
self.amp * (t_sec * self.omega).sin()
fn at(&self, t_sec: f32) -> Flt {
self.amp * (t_sec * self.omega).cast::<Flt>().sin()
}
}
impl TimeVarying3 for Sinusoid3 {
fn at(&self, t_sec: Flt) -> Vec3 {
fn at(&self, t_sec: f32) -> Vec3 {
self.amp * flt::Real::from_primitive(t_sec * self.omega).sin()
}
}
@@ -154,18 +153,18 @@ impl TimeVarying3 for Sinusoid3 {
#[derive(Clone)]
pub struct Gated<T> {
inner: T,
start: Flt,
end: Flt,
start: f32,
end: f32,
}
impl<T> Gated<T> {
pub fn new(inner: T, start: Flt, end: Flt) -> Self {
pub fn new(inner: T, start: f32, end: f32) -> Self {
Self { inner, start, end }
}
}
impl<T: TimeVarying1> TimeVarying1 for Gated<T> {
fn at(&self, t_sec: Flt) -> Flt {
fn at(&self, t_sec: f32) -> Flt {
if (self.start..self.end).contains(&t_sec) {
self.inner.at(t_sec)
} else {
@@ -175,7 +174,7 @@ impl<T: TimeVarying1> TimeVarying1 for Gated<T> {
}
impl<T: TimeVarying3> TimeVarying3 for Gated<T> {
fn at(&self, t_sec: Flt) -> Vec3 {
fn at(&self, t_sec: f32) -> Vec3 {
if (self.start..self.end).contains(&t_sec) {
self.inner.at(t_sec)
} else {
@@ -187,23 +186,23 @@ impl<T: TimeVarying3> TimeVarying3 for Gated<T> {
#[derive(Clone)]
pub struct Shifted<T> {
inner: T,
start_at: Flt,
start_at: f32,
}
impl<T> Shifted<T> {
pub fn new(inner: T, start_at: Flt) -> Self {
pub fn new(inner: T, start_at: f32) -> Self {
Self { inner, start_at }
}
}
impl<T: TimeVarying1> TimeVarying1 for Shifted<T> {
fn at(&self, t_sec: Flt) -> Flt {
fn at(&self, t_sec: f32) -> Flt {
self.inner.at(t_sec - self.start_at)
}
}
impl<T: TimeVarying3> TimeVarying3 for Shifted<T> {
fn at(&self, t_sec: Flt) -> Vec3 {
fn at(&self, t_sec: f32) -> Vec3 {
self.inner.at(t_sec - self.start_at)
}
}