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

Convert Stimulus to use f32 instead of Flt

Browse Source 
It's returned in a density format, so the scaling sensitive part
can be done internally by SimState. It should give pretty similar
results.
This commit is contained in:
Colin
2021-06-07 18:36:15 -07:00
parent 8779371f42
commit b97c298573
5 changed files with 42 additions and 43 deletions
Download Patch File Download Diff File Expand all files Collapse all files

4
examples/wrapped_torus.rs
View File

@@ -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);

10
src/bin/pml_tuning.rs
View File

@@ -27,20 +27,20 @@ struct PmlStim<F> {
t_end: f32,
feat_size: f32,
}
impl<F: Fn(Index) -> (Vec3, f32) + Sync> AbstractStimulus for PmlStim<F> {
fn at(&self, t_sec: f32, pos: Meters) -> Vec3 {
impl<F: Fn(Index) -> (Vec3<f32>, f32) + Sync> AbstractStimulus for PmlStim<F> {
fn at(&self, t_sec: f32, pos: Meters) -> Vec3<f32> {
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<F>(state: &mut StaticSim, frames: u32, f: F)
where F: Fn(Index) -> (Vec3, f32) + Sync // returns (E vector, omega)
where F: Fn(Index) -> (Vec3<f32>, 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<R, F>(state: &mut StaticSim, frames: u32, reg: R, f: F)
where
R: Region,
F: Fn(Index) -> (Vec3, f32) + Sync,
F: Fn(Index) -> (Vec3<f32>, f32) + Sync,
{
let feat = state.feature_size();
apply_stim_full_interior(state, frames, |idx| {

6
src/driver.rs
View File

@@ -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<f32> {
self.stim.at(t_sec, pos) * (self.frame_interval as f32)
}
}

28
src/sim.rs
View File

@@ -461,7 +461,7 @@ impl<M: Material + Clone + 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).cast::<Real>() * Real::from_primitive(timestep);
*e += value;
});
trace!("apply_stimulus end");
@@ -1641,20 +1641,20 @@ mod test {
feat_size: f32,
sqrt_vol: f32,
}
impl<F: Fn(Index) -> (Vec3, f32) + Sync> AbstractStimulus for PmlStim<F> {
fn at(&self, t_sec: f32, pos: Meters) -> Vec3 {
impl<F: Fn(Index) -> (Vec3<f32>, f32) + Sync> AbstractStimulus for PmlStim<F> {
fn at(&self, t_sec: f32, pos: Meters) -> Vec3<f32> {
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<F, S: GenericSim>(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>, f32) + Sync // returns (E vector, cycles (like Hz))
{
let stim = PmlStim {
f,
@@ -1677,7 +1677,7 @@ mod test {
fn pml_test_over_region<R, F, S>(state: &mut S, reg: R, f: F) -> f32
where
R: Region,
F: Fn(Index) -> (Vec3, f32) + Sync,
F: Fn(Index) -> (Vec3<f32>, f32) + Sync,
S: GenericSim
{
let feat = state.feature_size();
@@ -1689,13 +1689,13 @@ mod test {
}
})
}
fn pml_test_uniform_over_region<R: Region, S: GenericSim>(state: &mut S, reg: R, e: Vec3, hz: f32) -> f32 {
fn pml_test_uniform_over_region<R: Region, S: GenericSim>(state: &mut S, reg: R, e: Vec3<f32>, 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) -> f32 {
fn pml_test_at<S: GenericSim>(state: &mut S, e: Vec3<f32>, center: Index) -> f32 {
pml_test_full_interior(state, |idx| {
if idx == center {
(e, 1.0)
@@ -1705,14 +1705,14 @@ mod test {
})
}
fn pml_test<S: GenericSim>(state: &mut S, e: Vec3) {
fn pml_test<S: GenericSim>(state: &mut S, e: Vec3<f32>) {
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<P: GenericSim, B: GenericSim>(pml_state: &mut P, baseline_state: &mut B, e: Vec3) {
fn pml_test_against_baseline<P: GenericSim, B: GenericSim>(pml_state: &mut P, baseline_state: &mut B, e: Vec3<f32>) {
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<F: Fn(Vec3) -> Vec3>(
size: Index, e: Vec3, shuffle: F
fn pml_ineffective_mono_linear_test<F: Fn(Vec3<f32>) -> Vec3<f32>>(
size: Index, e: Vec3<f32>, 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);

37
src/stim.rs
View File

@@ -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<f32>;
}
// impl<T: AbstractStimulus> AbstractStimulus for &T {
@@ -14,13 +13,13 @@ pub trait AbstractStimulus: Sync {
// }
impl<T: AbstractStimulus> AbstractStimulus for Vec<T> {
fn at(&self, t_sec: f32, pos: Meters) -> Vec3 {
fn at(&self, t_sec: f32, pos: Meters) -> Vec3<f32> {
self.iter().map(|s| s.at(t_sec, pos)).sum()
}
}
impl AbstractStimulus for Box<dyn AbstractStimulus> {
fn at(&self, t_sec: f32, pos: Meters) -> Vec3 {
fn at(&self, t_sec: f32, pos: Meters) -> Vec3<f32> {
(**self).at(t_sec, pos)
}
}
@@ -41,7 +40,7 @@ impl<R, T> Stimulus<R, T> {
}
impl<R: Region + Sync, T: TimeVarying3 + Sync> AbstractStimulus for Stimulus<R, T> {
fn at(&self, t_sec: f32, pos: Meters) -> Vec3 {
fn at(&self, t_sec: f32, pos: Meters) -> Vec3<f32> {
if self.region.contains(pos) {
self.stim.at(t_sec)
} else {
@@ -67,11 +66,11 @@ impl<R, T> CurlStimulus<R, T> {
}
impl<R: Region + Sync, T: TimeVarying1 + Sync> AbstractStimulus for CurlStimulus<R, T> {
fn at(&self, t_sec: f32, pos: Meters) -> Vec3 {
fn at(&self, t_sec: f32, pos: Meters) -> Vec3<f32> {
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<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: f32) -> Flt;
fn at(&self, t_sec: f32) -> f32;
fn shifted(self, new_start: f32) -> Shifted<Self> {
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<f32>;
fn shifted(self, new_start: f32) -> Shifted<Self> {
Shifted::new(self, new_start)
}
@@ -109,8 +108,8 @@ pub struct Sinusoid<A> {
omega: f32,
}
pub type Sinusoid1 = Sinusoid<Flt>;
pub type Sinusoid3 = Sinusoid<Vec3>;
pub type Sinusoid1 = Sinusoid<f32>;
pub type Sinusoid3 = Sinusoid<Vec3<f32>>;
impl<A> Sinusoid<A> {
pub fn new(amp: A, freq: f32) -> Self {
@@ -139,14 +138,14 @@ impl<A> Sinusoid<A> {
}
impl TimeVarying1 for Sinusoid1 {
fn at(&self, t_sec: f32) -> Flt {
self.amp * (t_sec * self.omega).cast::<Flt>().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<f32> {
self.amp * (t_sec * self.omega).sin()
}
}
@@ -164,7 +163,7 @@ impl<T> Gated<T> {
}
impl<T: TimeVarying1> TimeVarying1 for Gated<T> {
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<T: TimeVarying1> TimeVarying1 for Gated<T> {
}
impl<T: TimeVarying3> TimeVarying3 for Gated<T> {
fn at(&self, t_sec: f32) -> Vec3 {
fn at(&self, t_sec: f32) -> Vec3<f32> {
if (self.start..self.end).contains(&t_sec) {
self.inner.at(t_sec)
} else {
@@ -196,13 +195,13 @@ impl<T> Shifted<T> {
}
impl<T: TimeVarying1> TimeVarying1 for Shifted<T> {
fn at(&self, t_sec: f32) -> Flt {
fn at(&self, t_sec: f32) -> f32 {
self.inner.at(t_sec - self.start_at)
}
}
impl<T: TimeVarying3> TimeVarying3 for Shifted<T> {
fn at(&self, t_sec: f32) -> Vec3 {
fn at(&self, t_sec: f32) -> Vec3<f32> {
self.inner.at(t_sec - self.start_at)
}
}