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

Stimulus: change at method to accept feat_size: f32, loc: Index

Browse Source
This commit is contained in:
colin
2022-08-18 16:21:21 -07:00
parent 0fa8cb0d20
commit cf2d21f975
4 changed files with 73 additions and 66 deletions
Download Patch File Download Diff File Expand all files Collapse all files

11
crates/coremem/src/geom/units.rs
View File

@@ -188,6 +188,17 @@ impl Index {
}
}
impl Into<Vec3u> for Index {
fn into(self) -> Vec3u {
self.0
}
}
impl From<Vec3u> for Index {
fn from(v: Vec3u) -> Self {
Self(v)
}
}
impl Coord for Index {
fn to_meters(&self, feature_size: f32) -> Meters {
Meters(Vec3::from(self.0) * feature_size)

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

@@ -304,8 +304,8 @@ impl<R: Real, M: Material<R> + cross::mat::Material<R> + Send + Sync> SimState<R
trace!("apply_stimulus begin");
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 densities = stim.at(t_sec, pos_meters);
let idx = Index::new(x as u32, y as u32, z as u32);
let densities = stim.at(t_sec, feature_size, idx);
let value_e = densities.e.cast::<R>() * timestep.cast::<R>();
let value_h = densities.h.cast::<R>() * timestep.cast::<R>();
*e += value_e;
@@ -792,7 +792,7 @@ impl<R: Real> CellState<R> {
#[cfg(test)]
mod test {
use super::*;
use crate::geom::{Cube, Meters, Region, WorldRegion};
use crate::geom::{Cube, Region, WorldRegion};
use crate::sim::legacy::mat::{AnisomorphicConductor, IsomorphicConductor, Pml, Static};
use crate::stim::Fields;
use crate::real::{R64, ToFloat as _};
@@ -1170,13 +1170,12 @@ mod test {
struct ConstStim {
loc: Index,
feature_size: f32,
e: Vec3<f32>,
h: Vec3<f32>,
}
impl Stimulus for ConstStim {
fn at(&self, _t_sec: f32, loc: Meters) -> Fields {
if loc.to_index(self.feature_size) == self.loc {
fn at(&self, _t_sec: f32, _feat_size: f32, loc: Index) -> Fields {
if loc == self.loc {
Fields { e: self.e, h: self.h }
} else {
Default::default()
@@ -1194,7 +1193,6 @@ mod test {
//println!("{}", g.exp());
state.step_multiple(1, &ConstStim {
loc: Index::new(0, 0, 1000),
feature_size: 1e-6,
e: Vec3::new_z(gate*sig/state.timestep()),
h: Vec3::zero(),
});
@@ -1213,7 +1211,6 @@ mod test {
let gate = 1e9 * 0.5*(1.0 - angle.cos());
state.step_multiple(1, &ConstStim {
loc: Index::new(10, 10, 10),
feature_size: 1e-6,
e: Vec3::new_z(gate*sig/state.timestep()),
h: Vec3::zero(),
});
@@ -1236,7 +1233,6 @@ mod test {
let gate = (progress*f32::pi()).sin();
state.step_multiple(1, &ConstStim {
loc: Index::new(100, 0, 0),
feature_size: 1e-6,
e: Vec3::new_z(gate*sig/state.timestep()),
h: Vec3::zero(),
});
@@ -1312,14 +1308,13 @@ mod test {
/// Maps location -> (field direction, frequency)
f: F,
t_end: f32,
feat_size: f32,
sqrt_vol: f32,
}
impl<F: Fn(Index) -> (Vec3<f32>, f32) + Sync> Stimulus for PmlStim<F> {
fn at(&self, t_sec: f32, pos: Meters) -> Fields {
fn at(&self, t_sec: f32, _feat_size: f32, loc: Index) -> 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 (e, hz) = (self.f)(loc);
let sig_angle = angle*hz;
let sig = sig_angle.sin();
// TODO: test H component?
@@ -1337,7 +1332,6 @@ mod test {
let stim = PmlStim {
f,
t_end: 50.0 * state.timestep(),
feat_size: state.feature_size(),
sqrt_vol: state.volume().sqrt(),
};
for _t in 0..50 {

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

@@ -2,7 +2,7 @@ use serde::{Deserialize, Serialize};
use log::{info, trace, warn};
use crate::diagnostics::SyncDiagnostics;
use crate::geom::{Coord, Index};
use crate::geom::Index;
use crate::real::Real;
use crate::sim::{AbstractSim, Fields, GenericSim, StaticSim};
use crate::stim::Stimulus;
@@ -269,6 +269,7 @@ where
}
}
// used for test
fn iterate_stim<R, S, F>(
stim: &S, dim: Index, feature_size: f32, t_sec: f32, timestep: R, mut f: F
)
@@ -282,8 +283,7 @@ where
for y in 0..dim.y() {
for x in 0..dim.x() {
let pos_idx = Index::new(x, y, z);
let pos_meters = pos_idx.to_meters(feature_size);
let densities = stim.at(t_sec, pos_meters);
let densities = stim.at(t_sec, feature_size, pos_idx);
let (value_e, value_h) = (densities.e.cast::<R>() * timestep, densities.h.cast::<R>() * timestep);
f(pos_idx, value_e, value_h)
}

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

@@ -3,6 +3,7 @@ use crate::cross::vec::Vec3;
use crate::geom::{Coord as _, HasCrossSection, Index, Meters, Region};
use coremem_cross::compound::list::{Visit, Visitor};
use coremem_cross::dim::OffsetDimSlice;
use coremem_cross::vec::Vec3u;
use rand;
/// field densities
@@ -70,32 +71,48 @@ impl std::ops::Mul<f32> for FieldMags {
}
pub trait Stimulus: Sync {
// 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;
fn at(&self, t_sec: f32, feat_size: f32, loc: Index) -> Fields {
let mut dat = OffsetDimSlice::new(loc.into(), Vec3u::new(1, 1, 1), [Fields::default()]);
self.eval_into(t_sec, feat_size, 1.0, dat.as_mut());
let [fields] = dat.into_inner();
fields
}
// TODO: could remove the `scale` param if we parameterized this array over some `F: FnMut(Fields)` instead `Fields`.
// that would also allow easier unzipping in the SpirvSim code.
/// bulk version of at. evaluate several positions at once and populate the output region.
fn eval_into(&self, t_sec: f32, feat_size: f32, scale: f32, into: OffsetDimSlice<&mut [Fields]>) {
for (idx, out) in into.enumerated() {
*out += self.at(t_sec, Index(idx).to_meters(feat_size)) * scale;
*out += self.at(t_sec, feat_size, Index::from(idx)) * scale;
}
}
}
/// used as a MapVisitor in order to evaluate each Stimulus in a List at a specific time/place.
struct StimulusEvaluator {
fields: Fields,
t_sec: f32,
pos: Meters,
}
impl<S: Stimulus> Visitor<&S> for &mut StimulusEvaluator {
fn visit(&mut self, next: &S) {
self.fields += next.at(self.t_sec, self.pos);
}
}
// struct StimulusEvaluator {
// fields: Fields,
// t_sec: f32,
// feat_size: f32,
// loc: Index,
// }
//
// impl<S: Stimulus> Visitor<&S> for &mut StimulusEvaluator {
// fn visit(&mut self, next: &S) {
// self.fields += next.at(self.t_sec, self.feat_size, self.loc);
// }
// }
//
// impl<L: Sync> Stimulus for L
// where
// for<'a, 'b> &'a L: Visit<&'b mut StimulusEvaluator>,
// {
// fn at(&self, t_sec: f32, pos: Meters) -> Fields {
// let mut ev = StimulusEvaluator { t_sec, pos, fields: Fields::default()};
// self.visit(&mut ev);
// ev.fields
// }
// }
struct StimEvalInto<'a> {
into: OffsetDimSlice<&'a mut [Fields]>,
@@ -110,26 +127,10 @@ impl<'a, S: Stimulus> Visitor<&S> for StimEvalInto<'a> {
}
}
// impl<L: Sync> Stimulus for L
// where
// for<'a, 'b> &'a L: Visit<&'b mut StimulusEvaluator>,
// {
// fn at(&self, t_sec: f32, pos: Meters) -> Fields {
// let mut ev = StimulusEvaluator { t_sec, pos, fields: Fields::default()};
// self.visit(&mut ev);
// ev.fields
// }
// }
impl<L: Sync> Stimulus for L
where
for<'a, 'b> &'a L: Visit<StimEvalInto<'b>>,
{
fn at(&self, _t_sec: f32, _pos: Meters) -> Fields {
// TODO: if we replace `_pos` with `feat_size`, `Index`, we can replace this with a generic
// impl that calls eval_into one a one-len slice.
unimplemented!();
}
fn eval_into(&self, t_sec: f32, feat_size: f32, scale: f32, into: OffsetDimSlice<&mut [Fields]>) {
let ev = StimEvalInto { t_sec, feat_size, scale, into };
self.visit(ev);
@@ -159,8 +160,8 @@ impl DynStimuli {
}
}
impl Stimulus for DynStimuli {
fn at(&self, t_sec: f32, pos: Meters) -> Fields {
self.0.iter().map(|i| i.at(t_sec, pos))
fn at(&self, t_sec: f32, feat_size: f32, loc: Index) -> Fields {
self.0.iter().map(|i| i.at(t_sec, feat_size, loc))
.fold(Fields::default(), core::ops::Add::add)
}
fn eval_into(&self, t_sec: f32, feat_size: f32, scale: f32, mut into: OffsetDimSlice<&mut [Fields]>) {
@@ -178,7 +179,7 @@ impl Stimulus for DynStimuli {
pub struct NoopStimulus;
impl Stimulus for NoopStimulus {
fn at(&self, _t_sec: f32, _pos: Meters) -> Fields {
fn at(&self, _t_sec: f32, _feat_size: f32, _loc: Index) -> Fields {
Fields::default()
}
}
@@ -197,7 +198,7 @@ impl UniformStimulus {
}
impl Stimulus for UniformStimulus {
fn at(&self, _t_sec: f32, _pos: Meters) -> Fields {
fn at(&self, _t_sec: f32, _feat_size: f32, _loc: Index) -> Fields {
self.fields
}
}
@@ -233,10 +234,10 @@ impl RngStimulus {
}
impl Stimulus for RngStimulus {
fn at(&self, t_sec: f32, pos: Meters) -> Fields {
fn at(&self, t_sec: f32, feat_size: f32, loc: Index) -> Fields {
Fields {
e: self.gen(t_sec, pos, self.e_scale, 0),
h: self.gen(t_sec, pos, self.h_scale, 0x7de3)
e: self.gen(t_sec, loc.to_meters(feat_size), self.e_scale, 0),
h: self.gen(t_sec, loc.to_meters(feat_size), self.h_scale, 0x7de3)
}
}
}
@@ -257,9 +258,9 @@ impl<R, S> RegionGated<R, S> {
}
impl<R: Region + Sync, S: Stimulus + Sync> Stimulus for RegionGated<R, S> {
fn at(&self, t_sec: f32, pos: Meters) -> Fields {
if self.region.contains(pos) {
self.stim.at(t_sec, pos)
fn at(&self, t_sec: f32, feat_size: f32, loc: Index) -> Fields {
if self.region.contains(loc.to_meters(feat_size)) {
self.stim.at(t_sec, feat_size, loc)
} else {
Fields::default()
}
@@ -281,7 +282,8 @@ impl<R, T> CurlStimulus<R, T> {
}
impl<R: Region + HasCrossSection, T: TimeVarying + Sync> Stimulus for CurlStimulus<R, T> {
fn at(&self, t_sec: f32, pos: Meters) -> Fields {
fn at(&self, t_sec: f32, feat_size: f32, loc: Index) -> Fields {
let pos = loc.to_meters(feat_size);
if self.region.contains(pos) {
let FieldMags { e, h } = self.stim.at(t_sec);
let rotational = self.region.cross_section_normal(pos);
@@ -395,9 +397,9 @@ impl<A: TimeVarying> TimeVarying for Exp<A> {
}
impl<A: Stimulus> Stimulus for Exp<A> {
fn at(&self, t_sec: f32, pos: Meters) -> Fields {
fn at(&self, t_sec: f32, feat_size: f32, loc: Index) -> Fields {
let scale = (t_sec * -self.tau).exp();
let inner = self.amp.at(t_sec, pos);
let inner = self.amp.at(t_sec, feat_size, loc);
Fields {
e: inner.e * scale,
h: inner.h * scale,
@@ -433,9 +435,9 @@ impl<T: TimeVarying> TimeVarying for Gated<T> {
}
impl<T: Stimulus> Stimulus for Gated<T> {
fn at(&self, t_sec: f32, pos: Meters) -> Fields {
fn at(&self, t_sec: f32, feat_size: f32, loc: Index) -> Fields {
if (self.start..self.end).contains(&t_sec) {
self.inner.at(t_sec, pos)
self.inner.at(t_sec, feat_size, loc)
} else {
Default::default()
}
@@ -466,8 +468,8 @@ impl<T: TimeVarying> TimeVarying for Shifted<T> {
}
impl<T: Stimulus> Stimulus for Shifted<T> {
fn at(&self, t_sec: f32, pos: Meters) -> Fields {
self.inner.at(t_sec - self.start_at, pos)
fn at(&self, t_sec: f32, feat_size: f32, loc: Index) -> Fields {
self.inner.at(t_sec - self.start_at, feat_size, loc)
}
fn eval_into(&self, t_sec: f32, feat_size: f32, scale: f32, into: OffsetDimSlice<&mut [Fields]>) {
self.inner.eval_into(t_sec - self.start_at, feat_size, scale, into)
@@ -529,7 +531,7 @@ mod test {
};
// stim acts as e: 1.0, h: 0.0
let stim = CurlStimulus::new(region, 1.0);
assert_eq!(stim.at(0.0, Meters::new(0.0, 0.0, 0.0)), Fields {
assert_eq!(stim.at(0.0, 1.0, Index::new(0, 0, 0)), Fields {
e: Vec3::new(1.0, 0.0, 0.0),
h: Vec3::new(0.0, 0.0, 0.0),
});
@@ -544,7 +546,7 @@ mod test {
};
// stim acts as e: -3.0, h: 0.0
let stim = CurlStimulus::new(region, -3.0);
assert_eq!(stim.at(0.0, Meters::new(0.0, 0.0, 0.0)), Fields {
assert_eq!(stim.at(0.0, 1.0, Index::new(0, 0, 0)), Fields {
e: Vec3::new(-3.0, 0.0, 0.0),
h: Vec3::new(0.0, 0.0, 0.0),
});
@@ -557,7 +559,7 @@ mod test {
};
// stim acts as e: 2.0, h: 0.0
let stim = CurlStimulus::new(region, 2.0);
let Fields { e, h: _ } = stim.at(0.0, Meters::new(0.0, 0.0, 0.0));
let Fields { e, h: _ } = stim.at(0.0, 1.0, Index::new(0, 0, 0));
assert!(e.distance(
Vec3::new(0.0, -1.0, 1.0).with_mag(2.0).unwrap()
) < 1e-6