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

stim: strongly-type the return type of AbstractSim::at with a `Fields` struct 

this will help me not mix up the E and H fields.
This commit is contained in:
colin 2022-07-30 17:17:17 -07:00
parent 6a511943f7
commit 4361167f99
4 changed files with 77 additions and 56 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

5
crates/coremem/src/driver.rs
View File

@ -5,8 +5,7 @@ use crate::real::Real;
use crate::render::{self, MultiRenderer, Renderer};
use crate::sim::AbstractSim;
use crate::sim::units::{Frame, Time};
use crate::stim::AbstractStimulus;
use crate::cross::vec::Vec3;
use crate::stim::{self, AbstractStimulus};
use log::{info, trace};
use serde::{Deserialize, Serialize};
@ -287,7 +286,7 @@ struct StimuliAdapter {
}
impl AbstractStimulus for StimuliAdapter {
fn at(&self, t_sec: f32, pos: Meters) -> (Vec3<f32>, Vec3<f32>) {
fn at(&self, t_sec: f32, pos: Meters) -> stim::Fields {
self.stim.at(t_sec, pos)
// TODO: remove this stuff (here only for testing)
/*

18
crates/coremem/src/sim/legacy/mod.rs
View File

@ -305,9 +305,9 @@ impl<R: Real, M: Material<R> + cross::mat::Material<R> + Send + Sync> SimState<R
Zip::from(ndarray::indices_of(&self.e)).and(&mut self.e).and(&mut self.h).par_for_each(
|(z, y, x), e, h| {
let pos_meters = Index((x as u32, y as u32, z as u32).into()).to_meters(feature_size);
let (density_e, density_h) = stim.at(t_sec, pos_meters);
let value_e = density_e.cast::<R>() * timestep.cast::<R>();
let value_h = density_h.cast::<R>() * timestep.cast::<R>();
let densities = stim.at(t_sec, pos_meters);
let value_e = densities.e.cast::<R>() * timestep.cast::<R>();
let value_h = densities.h.cast::<R>() * timestep.cast::<R>();
*e += value_e;
*h += value_h;
});
@ -794,6 +794,7 @@ mod test {
use super::*;
use crate::geom::{Cube, Meters, Region, WorldRegion};
use crate::sim::legacy::mat::{AnisomorphicConductor, IsomorphicConductor, Pml, Static};
use crate::stim::Fields;
use crate::real::{R64, ToFloat as _};
use float_eq::assert_float_eq;
use more_asserts::*;
@ -1173,9 +1174,9 @@ mod test {
h: Vec3<f32>,
}
impl AbstractStimulus for ConstStim {
fn at(&self, _t_sec: f32, loc: Meters) -> (Vec3<f32>, Vec3<f32>) {
fn at(&self, _t_sec: f32, loc: Meters) -> Fields {
if loc.to_index(self.feature_size) == self.loc {
(self.e, self.h)
Fields { e: self.e, h: self.h }
} else {
Default::default()
}
@ -1314,14 +1315,17 @@ mod test {
sqrt_vol: f32,
}
impl<F: Fn(Index) -> (Vec3<f32>, f32) + Sync> AbstractStimulus for PmlStim<F> {
fn at(&self, t_sec: f32, pos: Meters) -> (Vec3<f32>, Vec3<f32>) {
fn at(&self, t_sec: f32, pos: Meters) -> Fields {
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();
// TODO: test H component?
(e*(gate*sig / (self.t_end * self.sqrt_vol)), Vec3::zero())
Fields {
e: e*(gate*sig / (self.t_end * self.sqrt_vol)),
h: Vec3::zero()
}
}
}

12
crates/coremem/src/sim/spirv/mod.rs
View File

@ -228,16 +228,16 @@ where
)).par_map_collect(|(z, y, x)| {
let pos_idx = Index::new(x as _, y as _, z as _);
let pos_meters = pos_idx.to_meters(feature_size);
let (density_e, _density_h) = stim.at(t_sec, pos_meters);
density_e.cast::<R>() * timestep
let densities = stim.at(t_sec, pos_meters);
densities.e.cast::<R>() * timestep
});
let h = ndarray::Zip::from(ndarray::indices(
[dim.z() as usize, dim.y() as usize, dim.x() as usize]
)).par_map_collect(|(z, y, x)| {
let pos_idx = Index::new(x as _, y as _, z as _);
let pos_meters = pos_idx.to_meters(feature_size);
let (_density_e, density_h) = stim.at(t_sec, pos_meters);
density_h.cast::<R>() * timestep
let densities = stim.at(t_sec, pos_meters);
densities.h.cast::<R>() * timestep
});
trace!("eval_stimulus end");
(e, h)
@ -258,8 +258,8 @@ where
for x in 0..dim.x() {
let pos_idx = Index::new(x, y, z);
let pos_meters = pos_idx.to_meters(feature_size);
let (density_e, density_h) = stim.at(t_sec, pos_meters);
let (value_e, value_h) = (density_e.cast::<R>() * timestep, density_h.cast::<R>() * timestep);
let densities = stim.at(t_sec, pos_meters);
let (value_e, value_h) = (densities.e.cast::<R>() * timestep, densities.h.cast::<R>() * timestep);
f(pos_idx, value_e, value_h)
}
}

98
crates/coremem/src/stim.rs
View File

@ -3,10 +3,22 @@ use crate::cross::vec::Vec3;
use crate::geom::{Meters, Region};
use rand;
type Fields = (Vec3<f32>, Vec3<f32>);
/// field densities
#[derive(Clone, Copy, Debug, Default, PartialEq)]
pub struct Fields {
pub e: Vec3<f32>,
pub h: Vec3<f32>,
}
/// field magnitude densities (really, signed magnitude)
#[derive(Clone, Copy, Debug, Default, PartialEq)]
pub struct FieldMags {
pub e: f32,
pub h: f32,
}
pub trait AbstractStimulus: Sync {
// TODO: might be cleaner to return some `Fields` type instead of a tuple
// TODO: using floats for time and position is not great.
/// Return the (E, H) field which should be added PER-SECOND to the provided position/time.
fn at(&self, t_sec: f32, pos: Meters) -> Fields;
}
@ -19,13 +31,13 @@ pub trait AbstractStimulus: Sync {
impl<T: AbstractStimulus> AbstractStimulus for Vec<T> {
fn at(&self, t_sec: f32, pos: Meters) -> Fields {
let (mut e, mut h) = Fields::default();
let mut fields = Fields::default();
for i in self {
let (de, dh) = i.at(t_sec, pos);
e += de;
h += dh;
let Fields { e, h } = i.at(t_sec, pos);
fields.e += e;
fields.h += h;
}
(e, h)
fields
}
}
@ -64,7 +76,7 @@ impl AbstractStimulus for UniformStimulus {
impl TimeVarying for UniformStimulus {}
impl TimeVarying3 for UniformStimulus {
fn at(&self, _t_sec: f32) -> Fields {
(self.e, self.h)
Fields { e: self.e, h: self.h }
}
}
@ -100,7 +112,10 @@ impl RngStimulus {
impl AbstractStimulus for RngStimulus {
fn at(&self, t_sec: f32, pos: Meters) -> Fields {
(self.gen(t_sec, pos, self.e_scale, 0), self.gen(t_sec, pos, self.h_scale, 0x7de3))
Fields {
e: self.gen(t_sec, pos, self.e_scale, 0),
h: self.gen(t_sec, pos, self.h_scale, 0x7de3)
}
}
}
@ -148,13 +163,16 @@ 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) -> Fields {
if self.region.contains(pos) {
let (amt_e, amt_h) = self.stim.at(t_sec);
let FieldMags { e, h } = self.stim.at(t_sec);
let from_center_to_point = *pos - *self.center;
// TODO: is this inverted?
let rotational = from_center_to_point.cross(*self.axis);
let impulse_e = rotational.with_mag(amt_e.cast()).unwrap_or_default();
let impulse_h = rotational.with_mag(amt_h.cast()).unwrap_or_default();
(impulse_e, impulse_h)
let impulse_e = rotational.with_mag(e.cast()).unwrap_or_default();
let impulse_h = rotational.with_mag(h.cast()).unwrap_or_default();
Fields {
e: impulse_e,
h: impulse_h,
}
} else {
Fields::default()
}
@ -172,7 +190,7 @@ pub trait TimeVarying: Sized {
pub trait TimeVarying1: TimeVarying {
/// Retrieve the (E, H) impulse to apply PER-SECOND at the provided time (in seconds).
fn at(&self, t_sec: f32) -> (f32, f32);
fn at(&self, t_sec: f32) -> FieldMags;
}
pub trait TimeVarying3: TimeVarying {
@ -183,8 +201,8 @@ pub trait TimeVarying3: TimeVarying {
// assumed to represent the E field
impl TimeVarying for f32 {}
impl TimeVarying1 for f32 {
fn at(&self, _t_sec: f32) -> (f32, f32) {
(*self, 0.0)
fn at(&self, _t_sec: f32) -> FieldMags {
FieldMags { e: *self, h: 0.0 }
}
}
@ -227,20 +245,20 @@ impl<A> Sinusoid<A> {
impl<A> TimeVarying for Sinusoid<A> {}
impl TimeVarying1 for Sinusoid1 {
fn at(&self, t_sec: f32) -> (f32, f32) {
(
self.amp * (t_sec * self.omega).sin(),
0.0,
)
fn at(&self, t_sec: f32) -> FieldMags {
FieldMags {
e: self.amp * (t_sec * self.omega).sin(),
h: 0.0,
}
}
}
impl TimeVarying3 for Sinusoid3 {
fn at(&self, t_sec: f32) -> Fields {
(
self.amp * (t_sec * self.omega).sin(),
Vec3::zero(),
)
Fields {
e: self.amp * (t_sec * self.omega).sin(),
h: Vec3::zero(),
}
}
}
@ -269,20 +287,20 @@ impl<A> Exp<A> {
impl<A> TimeVarying for Exp<A> {}
impl TimeVarying1 for Exp1 {
fn at(&self, t_sec: f32) -> (f32, f32) {
(
self.amp * (t_sec * -self.tau).exp(),
0.0,
)
fn at(&self, t_sec: f32) -> FieldMags {
FieldMags {
e: self.amp * (t_sec * -self.tau).exp(),
h: 0.0,
}
}
}
impl TimeVarying3 for Exp3 {
fn at(&self, t_sec: f32) -> Fields {
(
self.amp * (t_sec * -self.tau).exp(),
Vec3::zero(),
)
Fields {
e: self.amp * (t_sec * -self.tau).exp(),
h: Vec3::zero(),
}
}
}
@ -301,7 +319,7 @@ impl<T> Gated<T> {
impl<T> TimeVarying for Gated<T> {}
impl<T: TimeVarying1> TimeVarying1 for Gated<T> {
fn at(&self, t_sec: f32) -> (f32, f32) {
fn at(&self, t_sec: f32) -> FieldMags {
if (self.start..self.end).contains(&t_sec) {
self.inner.at(t_sec)
} else {
@ -334,7 +352,7 @@ impl<T> Shifted<T> {
impl<T> TimeVarying for Shifted<T> {}
impl<T: TimeVarying1> TimeVarying1 for Shifted<T> {
fn at(&self, t_sec: f32) -> (f32, f32) {
fn at(&self, t_sec: f32) -> FieldMags {
self.inner.at(t_sec - self.start_at)
}
}
@ -354,9 +372,9 @@ mod test {
let x = $x;
let e = $e;
let h = $h;
let diff_e = (x.0 - e).mag();
let diff_e = (x.e - e).mag();
assert!(diff_e <= 0.001, "{:?} != {:?}", x, e);
let diff_h = (x.1 - h).mag();
let diff_h = (x.h - h).mag();
assert!(diff_h <= 0.001, "{:?} != {:?}", x, h);
}
}
@ -364,7 +382,7 @@ mod test {
#[test]
fn sinusoid3() {
let s = Sinusoid3::new(Vec3::new(10.0, 1.0, -100.0), 1000.0);
assert_eq!(s.at(0.0), (Vec3::zero(), Vec3::zero()));
assert_eq!(s.at(0.0), Fields::default());
assert_approx_eq!(s.at(0.00025),
Vec3::new(10.0, 1.0, -100.0), Vec3::zero()
);
@ -375,7 +393,7 @@ mod test {
#[test]
fn sinusoid3_from_wavelength() {
let s = Sinusoid3::from_wavelength(Vec3::new(10.0, 1.0, -100.0), 0.001);
assert_eq!(s.at(0.0), (Vec3::zero(), Vec3::zero()));
assert_eq!(s.at(0.0), Fields::default());
assert_approx_eq!(s.at(0.00025), Vec3::new(10.0, 1.0, -100.0), Vec3::zero());
assert_approx_eq!(s.at(0.00050), Vec3::zero(), Vec3::zero());
assert_approx_eq!(s.at(0.00075), Vec3::new(-10.0, -1.0, 100.0), Vec3::zero());