From b97c298573ea74e6be1aa8f235d86d1dfc93a6f4 Mon Sep 17 00:00:00 2001 From: Colin Date: Mon, 7 Jun 2021 18:36:15 -0700 Subject: [PATCH] Convert Stimulus to use f32 instead of Flt It's returned in a density format, so the scaling sensitive part can be done internally by SimState. It should give pretty similar results. --- examples/wrapped_torus.rs | 4 ++-- src/bin/pml_tuning.rs | 10 +++++----- src/driver.rs | 6 +++--- src/sim.rs | 28 ++++++++++++++-------------- src/stim.rs | 37 ++++++++++++++++++------------------- 5 files changed, 42 insertions(+), 43 deletions(-) diff --git a/examples/wrapped_torus.rs b/examples/wrapped_torus.rs index 4f806cc..a7d5afc 100644 --- a/examples/wrapped_torus.rs +++ b/examples/wrapped_torus.rs @@ -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 as Flt, current_duration * 2.0) + let pos_wave = Sinusoid1::from_wavelength(peak_stim, current_duration * 2.0) .half_cycle(); - let neg_wave = Sinusoid1::from_wavelength(-4.0*peak_stim as Flt, current_duration * 2.0) + let neg_wave = Sinusoid1::from_wavelength(-4.0*peak_stim, current_duration * 2.0) .half_cycle() .shifted(current_duration + current_break); diff --git a/src/bin/pml_tuning.rs b/src/bin/pml_tuning.rs index bd02d6d..768342d 100644 --- a/src/bin/pml_tuning.rs +++ b/src/bin/pml_tuning.rs @@ -27,20 +27,20 @@ struct PmlStim { t_end: f32, feat_size: f32, } -impl (Vec3, f32) + Sync> AbstractStimulus for PmlStim { - fn at(&self, t_sec: f32, pos: Meters) -> Vec3 { +impl (Vec3, f32) + Sync> AbstractStimulus for PmlStim { + 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; let sig = sig_angle.sin(); - e*Real::from_primitive(gate*sig) + e*gate*sig } } /// Apply some stimulus, and then let it decay and measure the ratio of energy left in the system fn apply_stim_full_interior(state: &mut StaticSim, frames: u32, f: F) -where F: Fn(Index) -> (Vec3, f32) + Sync // returns (E vector, omega) +where F: Fn(Index) -> (Vec3, f32) + Sync // returns (E vector, omega) { let stim = PmlStim { f, @@ -56,7 +56,7 @@ where F: Fn(Index) -> (Vec3, f32) + Sync // returns (E vector, omega) fn apply_stim_over_region(state: &mut StaticSim, frames: u32, reg: R, f: F) where R: Region, - F: Fn(Index) -> (Vec3, f32) + Sync, + F: Fn(Index) -> (Vec3, f32) + Sync, { let feat = state.feature_size(); apply_stim_full_interior(state, frames, |idx| { diff --git a/src/driver.rs b/src/driver.rs index 20ef19b..969b203 100644 --- a/src/driver.rs +++ b/src/driver.rs @@ -1,4 +1,4 @@ -use crate::{flt, flt::Flt, mat}; +use crate::{flt::Flt, mat}; use crate::geom::{Coord, Meters, Region, Vec3}; use crate::mat::{GenericMaterial, Material, Pml}; use crate::meas::{self, AbstractMeasurement}; @@ -223,8 +223,8 @@ struct StimuliAdapter { } impl AbstractStimulus for StimuliAdapter { - fn at(&self, t_sec: f32, pos: Meters) -> Vec3 { - self.stim.at(t_sec, pos) * flt::Real::from_primitive(self.frame_interval) + fn at(&self, t_sec: f32, pos: Meters) -> Vec3 { + self.stim.at(t_sec, pos) * (self.frame_interval as f32) } } diff --git a/src/sim.rs b/src/sim.rs index 2dec8ec..5bb9a1a 100644 --- a/src/sim.rs +++ b/src/sim.rs @@ -461,7 +461,7 @@ impl SimState { 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).cast::() * Real::from_primitive(timestep); *e += value; }); trace!("apply_stimulus end"); @@ -1641,20 +1641,20 @@ mod test { feat_size: f32, sqrt_vol: f32, } - impl (Vec3, f32) + Sync> AbstractStimulus for PmlStim { - fn at(&self, t_sec: f32, pos: Meters) -> Vec3 { + impl (Vec3, f32) + Sync> AbstractStimulus for PmlStim { + 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; let sig = sig_angle.sin(); - e*Real::from_primitive(gate*sig / (self.t_end * self.sqrt_vol)) + e*(gate*sig / (self.t_end * self.sqrt_vol)) } } /// Returns the energy ratio (Eend / Estart) fn pml_test_full_interior(state: &mut S, f: F) -> f32 - where F: Fn(Index) -> (Vec3, f32) + Sync // returns (E vector, cycles (like Hz)) + where F: Fn(Index) -> (Vec3, f32) + Sync // returns (E vector, cycles (like Hz)) { let stim = PmlStim { f, @@ -1677,7 +1677,7 @@ mod test { fn pml_test_over_region(state: &mut S, reg: R, f: F) -> f32 where R: Region, - F: Fn(Index) -> (Vec3, f32) + Sync, + F: Fn(Index) -> (Vec3, f32) + Sync, S: GenericSim { let feat = state.feature_size(); @@ -1689,13 +1689,13 @@ mod test { } }) } - fn pml_test_uniform_over_region(state: &mut S, reg: R, e: Vec3, hz: f32) -> f32 { + fn pml_test_uniform_over_region(state: &mut S, reg: R, e: Vec3, hz: f32) -> f32 { pml_test_over_region(state, reg, |_idx| { (e, hz) }) } - fn pml_test_at(state: &mut S, e: Vec3, center: Index) -> f32 { + fn pml_test_at(state: &mut S, e: Vec3, center: Index) -> f32 { pml_test_full_interior(state, |idx| { if idx == center { (e, 1.0) @@ -1705,14 +1705,14 @@ mod test { }) } - fn pml_test(state: &mut S, e: Vec3) { + fn pml_test(state: &mut S, e: Vec3) { let center = Index(state.size().0 / 2); let pml = pml_test_at(state, e, center); // PML should absorb all energy assert_float_eq!(pml, 0.0, abs <= 2e-5); } - fn pml_test_against_baseline(pml_state: &mut P, baseline_state: &mut B, e: Vec3) { + fn pml_test_against_baseline(pml_state: &mut P, baseline_state: &mut B, e: Vec3) { assert_eq!(pml_state.size(), baseline_state.size()); let center = Index(pml_state.size().0 / 2); let en_pml = pml_test_at(pml_state, e, center); @@ -1787,8 +1787,8 @@ mod test { } /// Despite PML existing, verify that it doesn't interfere with a specific wave. - fn pml_ineffective_mono_linear_test Vec3>( - size: Index, e: Vec3, shuffle: F + fn pml_ineffective_mono_linear_test) -> Vec3>( + size: Index, e: Vec3, shuffle: F ) { let mut state = SimState::new(size, 1e-6); let timestep = state.timestep(); @@ -1797,8 +1797,8 @@ mod test { 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.cast()); - Pml::new(coord_stretch) + let coord_stretch = shuffle(cs); + Pml::new(coord_stretch.cast()) }); let center = Index(state.size().0 / 2); let en = pml_test_at(&mut state, e, center); diff --git a/src/stim.rs b/src/stim.rs index 73b4cd9..966f366 100644 --- a/src/stim.rs +++ b/src/stim.rs @@ -1,10 +1,9 @@ -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: f32, pos: Meters) -> Vec3; + fn at(&self, t_sec: f32, pos: Meters) -> Vec3; } // impl AbstractStimulus for &T { @@ -14,13 +13,13 @@ pub trait AbstractStimulus: Sync { // } impl AbstractStimulus for Vec { - fn at(&self, t_sec: f32, 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 { - fn at(&self, t_sec: f32, pos: Meters) -> Vec3 { + fn at(&self, t_sec: f32, pos: Meters) -> Vec3 { (**self).at(t_sec, pos) } } @@ -41,7 +40,7 @@ impl Stimulus { } impl AbstractStimulus for Stimulus { - fn at(&self, t_sec: f32, pos: Meters) -> Vec3 { + fn at(&self, t_sec: f32, pos: Meters) -> Vec3 { if self.region.contains(pos) { self.stim.at(t_sec) } else { @@ -67,11 +66,11 @@ impl CurlStimulus { } impl AbstractStimulus for CurlStimulus { - fn at(&self, t_sec: f32, 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: Vec3 = from_center_to_point.cross(*self.axis).cast(); + let rotational = from_center_to_point.cross(*self.axis); let impulse = rotational.with_mag(amount.cast()); impulse } else { @@ -83,7 +82,7 @@ impl 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: f32) -> Flt; + fn at(&self, t_sec: f32) -> f32; fn shifted(self, new_start: f32) -> Shifted { Shifted::new(self, new_start) } @@ -94,7 +93,7 @@ pub trait TimeVarying1: Sized { pub trait TimeVarying3: Sized { /// Retrieve the E impulse to apply PER-SECOND at the provided time (in seconds). - fn at(&self, t_sec: f32) -> Vec3; + fn at(&self, t_sec: f32) -> Vec3; fn shifted(self, new_start: f32) -> Shifted { Shifted::new(self, new_start) } @@ -109,8 +108,8 @@ pub struct Sinusoid { omega: f32, } -pub type Sinusoid1 = Sinusoid; -pub type Sinusoid3 = Sinusoid; +pub type Sinusoid1 = Sinusoid; +pub type Sinusoid3 = Sinusoid>; impl Sinusoid { pub fn new(amp: A, freq: f32) -> Self { @@ -139,14 +138,14 @@ impl Sinusoid { } impl TimeVarying1 for Sinusoid1 { - fn at(&self, t_sec: f32) -> Flt { - self.amp * (t_sec * self.omega).cast::().sin() + fn at(&self, t_sec: f32) -> f32 { + self.amp * (t_sec * self.omega).sin() } } impl TimeVarying3 for Sinusoid3 { - fn at(&self, t_sec: f32) -> Vec3 { - self.amp * flt::Real::from_primitive(t_sec * self.omega).sin() + fn at(&self, t_sec: f32) -> Vec3 { + self.amp * (t_sec * self.omega).sin() } } @@ -164,7 +163,7 @@ impl Gated { } impl TimeVarying1 for Gated { - fn at(&self, t_sec: f32) -> Flt { + fn at(&self, t_sec: f32) -> f32 { if (self.start..self.end).contains(&t_sec) { self.inner.at(t_sec) } else { @@ -174,7 +173,7 @@ impl TimeVarying1 for Gated { } impl TimeVarying3 for Gated { - fn at(&self, t_sec: f32) -> Vec3 { + fn at(&self, t_sec: f32) -> Vec3 { if (self.start..self.end).contains(&t_sec) { self.inner.at(t_sec) } else { @@ -196,13 +195,13 @@ impl Shifted { } impl TimeVarying1 for Shifted { - fn at(&self, t_sec: f32) -> Flt { + fn at(&self, t_sec: f32) -> f32 { self.inner.at(t_sec - self.start_at) } } impl TimeVarying3 for Shifted { - fn at(&self, t_sec: f32) -> Vec3 { + fn at(&self, t_sec: f32) -> Vec3 { self.inner.at(t_sec - self.start_at) } }