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

Remove Real from src/sim.rs

Browse Source
This commit is contained in:
Colin
2021-06-09 15:28:22 -07:00
parent dfdbb43180
commit 490afbc4da
5 changed files with 51 additions and 47 deletions
Download Patch File Download Diff File Expand all files Collapse all files

1
src/mat.rs
View File

@@ -1,5 +1,4 @@
use crate::CellState; use crate::CellState;
use crate::flt::{self, Flt};
use crate::geom::{Line2d, Vec2, Vec3, Polygon2d}; use crate::geom::{Line2d, Vec2, Vec3, Polygon2d};
use crate::real::Real; use crate::real::Real;
use crate::sim::{PmlParameters, PmlState, StepParameters, StepParametersMut}; use crate::sim::{PmlParameters, PmlState, StepParameters, StepParametersMut};

8
src/post.rs
View File

@@ -1,5 +1,5 @@
//! Post-processing tools //! Post-processing tools
use crate::mat::GenericMaterial; use crate::mat::{GenericMaterial, Static};
use crate::meas::AbstractMeasurement; use crate::meas::AbstractMeasurement;
use crate::render::{ColorTermRenderer, Renderer as _, RenderConfig, SerializedFrame}; use crate::render::{ColorTermRenderer, Renderer as _, RenderConfig, SerializedFrame};
use crate::sim::{SimState, StaticSim}; use crate::sim::{SimState, StaticSim};
@@ -34,7 +34,7 @@ pub type Result<T> = std::result::Result<T, Error>;
pub struct Frame { pub struct Frame {
path: PathBuf, path: PathBuf,
data: SerializedFrame, data: SerializedFrame<StaticSim>,
} }
impl Frame { impl Frame {
@@ -111,9 +111,9 @@ impl Loader {
// decode to a valid but incorrect state... // decode to a valid but incorrect state...
let data = bincode::deserialize_from(&mut reader).or_else(|_| -> Result<_> { let data = bincode::deserialize_from(&mut reader).or_else(|_| -> Result<_> {
reader.seek(SeekFrom::Start(0)).unwrap(); reader.seek(SeekFrom::Start(0)).unwrap();
let data: SerializedFrame<SimState<GenericMaterial<crate::flt::Real>>> = let data: SerializedFrame<SimState<Static<f32>, f32>> =
bincode::deserialize_from(reader)?; bincode::deserialize_from(reader)?;
Ok(data.to_static()) Ok(SerializedFrame::to_static(data))
})?; })?;
Ok(Frame { Ok(Frame {
path: path.into(), path: path.into(),

10
src/real.rs
View File

@@ -4,6 +4,16 @@ use decorum::cmp::IntrinsicOrd;
pub use decorum::Real as decorum_Real; pub use decorum::Real as decorum_Real;
pub use num::Zero; pub use num::Zero;
// finite. non-nan, non-inf
pub mod checked {
pub use decorum::{R32, R64};
}
pub mod unchecked {
pub type R32 = f32;
pub type R64 = f64;
}
pub trait ToFloat { pub trait ToFloat {
fn to_f32(&self) -> f32 { fn to_f32(&self) -> f32 {
self.to_f64() as _ self.to_f64() as _

2
src/render.rs
View File

@@ -643,7 +643,7 @@ pub struct SerializedFrame<S=StaticSim> {
impl<S: GenericSim> SerializedFrame<S> { impl<S: GenericSim> SerializedFrame<S> {
pub fn to_static(self) -> SerializedFrame<StaticSim> { pub fn to_static(self) -> SerializedFrame<StaticSim> {
SerializedFrame { SerializedFrame {
state: self.state.to_static(), state: GenericSim::to_static(&self.state),
measurements: self.measurements, measurements: self.measurements,
} }
} }

77
src/sim.rs
View File

@@ -1,9 +1,7 @@
use crate::flt::Real;
use crate::geom::{Coord, Cube, Index, InvertedRegion, Meters, Region, Vec3, Vec3u}; use crate::geom::{Coord, Cube, Index, InvertedRegion, Meters, Region, Vec3, Vec3u};
use crate::mat::{self, GenericMaterial, Material, MaterialExt as _}; use crate::mat::{self, GenericMaterial, Material, MaterialExt as _};
use crate::real::{self, Real as _, ToFloat as _}; use crate::real::{checked, Real, ToFloat as _};
use crate::stim::AbstractStimulus; use crate::stim::AbstractStimulus;
use dyn_clone::{self, DynClone};
use log::trace; use log::trace;
use ndarray::{Array3, Zip}; use ndarray::{Array3, Zip};
use rayon::prelude::*; use rayon::prelude::*;
@@ -11,7 +9,7 @@ use serde::{Serialize, Deserialize};
use std::convert::From; use std::convert::From;
use std::iter::Sum; use std::iter::Sum;
pub type StaticSim = SimState<mat::Static<Real>, Real>; pub type StaticSim = SimState<mat::Static<f32>, f32>;
#[derive(Default, Copy, Clone, PartialEq, Serialize, Deserialize)] #[derive(Default, Copy, Clone, PartialEq, Serialize, Deserialize)]
pub struct PmlState<R> { pub struct PmlState<R> {
@@ -46,7 +44,7 @@ struct PmlAux<R> {
#[derive(Default, Copy, Clone, PartialEq, Serialize, Deserialize)] #[derive(Default, Copy, Clone, PartialEq, Serialize, Deserialize)]
struct ConvFlt<R>(R); struct ConvFlt<R>(R);
impl<R: real::Real> ConvFlt<R> { impl<R: Real> ConvFlt<R> {
///```tex ///```tex
/// Consider $v(t) = (a \delta(t) + b exp[-c t] u(t)) \ast y(t)$ /// Consider $v(t) = (a \delta(t) + b exp[-c t] u(t)) \ast y(t)$
/// This method will advance from self = $v(t)$ to self = $v(t + \Delta t)$ /// This method will advance from self = $v(t)$ to self = $v(t + \Delta t)$
@@ -114,8 +112,8 @@ pub struct StepParametersMut<'a, R> {
pml: Option<(&'a mut PmlState<R>, PmlParameters<R>)>, pml: Option<(&'a mut PmlState<R>, PmlParameters<R>)>,
} }
impl<'a, R: real::Real> StepParametersMut<'a, R> { impl<'a, R: Real> StepParametersMut<'a, R> {
pub fn new<R2: real::Real>( pub fn new<R2: Real>(
conductivity: Vec3<R2>, pml: Option<(&'a mut PmlState<R>, PmlParameters<R>)> conductivity: Vec3<R2>, pml: Option<(&'a mut PmlState<R>, PmlParameters<R>)>
) -> Self { ) -> Self {
Self { Self {
@@ -123,7 +121,7 @@ impl<'a, R: real::Real> StepParametersMut<'a, R> {
pml, pml,
} }
} }
pub fn with_conductivity<R2: real::Real>(mut self, conductivity: Vec3<R2>) -> Self { pub fn with_conductivity<R2: Real>(mut self, conductivity: Vec3<R2>) -> Self {
self.conductivity = conductivity.cast(); self.conductivity = conductivity.cast();
self self
} }
@@ -167,8 +165,8 @@ pub struct PmlParameters<R> {
pseudo_conductivity: Vec3<R>, pseudo_conductivity: Vec3<R>,
} }
impl<R: real::Real> PmlParameters<R> { impl<R: Real> PmlParameters<R> {
pub fn new<R2: real::Real>(pseudo_conductivity: Vec3<R2>) -> Self { pub fn new<R2: Real>(pseudo_conductivity: Vec3<R2>) -> Self {
Self { Self {
pseudo_conductivity: pseudo_conductivity.cast(), pseudo_conductivity: pseudo_conductivity.cast(),
} }
@@ -182,7 +180,7 @@ pub struct StepParameters<'a, R> {
pml: Option<(&'a PmlState<R>, PmlParameters<R>)>, pml: Option<(&'a PmlState<R>, PmlParameters<R>)>,
} }
impl<'a, R: real::Real> StepParameters<'a, R> { impl<'a, R: Real> StepParameters<'a, R> {
pub fn conductivity(&self) -> Vec3<R> { pub fn conductivity(&self) -> Vec3<R> {
self.conductivity self.conductivity
} }
@@ -200,7 +198,7 @@ impl<'a, R> From<StepParametersMut<'a, R>> for StepParameters<'a, R> {
} }
} }
pub trait GenericSim: Send + Sync + DynClone { pub trait GenericSim: Send + Sync {
fn sample(&self, pos: Meters) -> Sample; fn sample(&self, pos: Meters) -> Sample;
/// DEPRECATED. Use stimulus instead /// DEPRECATED. Use stimulus instead
fn impulse_e_meters(&mut self, pos: Meters, amount: Vec3<f32>); fn impulse_e_meters(&mut self, pos: Meters, amount: Vec3<f32>);
@@ -257,7 +255,6 @@ pub trait GenericSim: Send + Sync + DynClone {
state state
} }
} }
dyn_clone::clone_trait_object!(GenericSim);
impl<'a> dyn GenericSim + 'a { impl<'a> dyn GenericSim + 'a {
pub fn get<C: Coord>(&self, at: C) -> Sample { pub fn get<C: Coord>(&self, at: C) -> Sample {
@@ -372,7 +369,7 @@ impl<'a> dyn GenericSim + 'a {
} }
#[derive(Default, Clone, Serialize, Deserialize)] #[derive(Default, Clone, Serialize, Deserialize)]
pub struct SimState<M=GenericMaterial<Real>, R=Real> { pub struct SimState<M=GenericMaterial<checked::R32>, R=checked::R32> {
cells: Array3<M>, cells: Array3<M>,
e: Array3<Vec3<R>>, e: Array3<Vec3<R>>,
h: Array3<Vec3<R>>, h: Array3<Vec3<R>>,
@@ -400,7 +397,7 @@ impl<M, R> SimState<M, R> {
} }
} }
impl<M, R: real::Real> SimState<M, R> { impl<M, R: Real> SimState<M, R> {
pub fn feature_size(&self) -> f32 { pub fn feature_size(&self) -> f32 {
self.feature_size.to_f32() self.feature_size.to_f32()
} }
@@ -409,7 +406,7 @@ impl<M, R: real::Real> SimState<M, R> {
} }
} }
impl<R: real::Real, M: Default> SimState<M, R> { impl<R: Real, M: Default> SimState<M, R> {
pub fn new(size: Index, feature_size: f32) -> 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. // 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 // For 3d, we need 3 full timesteps in order to communicate information across the corner
@@ -431,7 +428,7 @@ impl<R: real::Real, M: Default> SimState<M, R> {
} }
} }
impl<R: real::Real, M: Material<R> + Clone + Send + Sync + 'static> SimState<M, R> { impl<R: Real, M: Material<R> + Send + Sync> SimState<M, R> {
pub fn step(&mut self) { pub fn step(&mut self) {
let time_step = self.timestep; let time_step = self.timestep;
let half_time_step = R::half() * time_step; let half_time_step = R::half() * time_step;
@@ -479,13 +476,13 @@ impl<R: real::Real, M: Material<R> + Clone + Send + Sync + 'static> SimState<M,
} }
} }
impl<R: real::Real, M: Material<R> + Clone> SimState<M, R> { impl<R: Real, M: Material<R> + Clone> SimState<M, R> {
pub fn fill_region<Reg: Region>(&mut self, region: &Reg, mat: M) { pub fn fill_region<Reg: Region>(&mut self, region: &Reg, mat: M) {
self.fill_region_using(region, |_idx: Index| mat.clone()); self.fill_region_using(region, |_idx: Index| mat.clone());
} }
} }
impl<R: real::Real, M: Material<R>> SimState<M, R> { impl<R: Real, M: Material<R>> SimState<M, R> {
pub fn fill_region_using<C, Reg, F>(&mut self, region: &Reg, f: F) pub fn fill_region_using<C, Reg, F>(&mut self, region: &Reg, f: F)
where where
Reg: Region, Reg: Region,
@@ -554,7 +551,7 @@ impl<R: real::Real, M: Material<R>> SimState<M, R> {
} }
} }
impl<R: real::Real, M: Material<R> + Clone + Send + Sync + 'static> GenericSim for SimState<M, R> { impl<R: Real, M: Material<R> + Send + Sync> GenericSim for SimState<M, R> {
fn sample(&self, pos: Meters) -> Sample { fn sample(&self, pos: Meters) -> Sample {
// TODO: smarter sampling than nearest neighbor? // TODO: smarter sampling than nearest neighbor?
let pos_sim = pos.to_index(self.feature_size()); let pos_sim = pos.to_index(self.feature_size());
@@ -610,7 +607,7 @@ impl<R: real::Real, M: Material<R> + Clone + Send + Sync + 'static> GenericSim f
} }
#[allow(unused)] #[allow(unused)]
impl<R: real::Real, M> SimState<M, R> { impl<R: Real, M> SimState<M, R> {
fn get_mat<C: Coord>(&self, c: C) -> &M { fn get_mat<C: Coord>(&self, c: C) -> &M {
let at = c.to_index(self.feature_size.cast()); let at = c.to_index(self.feature_size.cast());
&self.cells[at.row_major_idx()] &self.cells[at.row_major_idx()]
@@ -800,7 +797,7 @@ impl Sample {
/// Solve df/dt in: /// Solve df/dt in:
/// nabla_g = df_dt_coeff * df/dt + f_coeff * f + f_int_coeff * \int f + f_int_int_coeff * \int \int f /// nabla_g = df_dt_coeff * df/dt + f_coeff * f + f_int_coeff * \int f + f_int_int_coeff * \int \int f
#[allow(unused)] #[allow(unused)]
fn solve_step_diff_eq<R: real::Real>( fn solve_step_diff_eq<R: Real>(
nabla_g: Vec3<R>, nabla_g: Vec3<R>,
df_dt_coeff: Vec3<R>, df_dt_coeff: Vec3<R>,
f_coeff: Vec3<R>, f_coeff: Vec3<R>,
@@ -840,7 +837,7 @@ fn solve_step_diff_eq<R: real::Real>(
/// Solve df/dt in: /// Solve df/dt in:
/// nabla_g = df_dt_coeff * df/dt + f_coeff * f /// nabla_g = df_dt_coeff * df/dt + f_coeff * f
/// Lower-order (more efficient) version of `solve_step_diff_eq`. /// Lower-order (more efficient) version of `solve_step_diff_eq`.
fn solve_step_diff_eq_linear<R: real::Real>( fn solve_step_diff_eq_linear<R: Real>(
nabla_g: Vec3<R>, nabla_g: Vec3<R>,
df_dt_coeff: Vec3<R>, df_dt_coeff: Vec3<R>,
f_coeff: Vec3<R>, f_coeff: Vec3<R>,
@@ -868,7 +865,7 @@ struct ECell<'a, M, R> {
h: &'a Vec3<R>, h: &'a Vec3<R>,
} }
impl<'a, R: real::Real, M: Material<R>> HCell<'a, M, R> { impl<'a, R: Real, M: Material<R>> HCell<'a, M, R> {
pub fn b(&self) -> Vec3<R> { pub fn b(&self) -> Vec3<R> {
(*self.h + self.mat.m()) * R::mu0() (*self.h + self.mat.m()) * R::mu0()
} }
@@ -1039,7 +1036,7 @@ impl<'a, R: real::Real, M: Material<R>> HCell<'a, M, R> {
} }
} }
impl<'a, R: real::Real, M: Material<R>> ECell<'a, M, R> { impl<'a, R: Real, M: Material<R>> ECell<'a, M, R> {
/// delta_t = timestep covered by this function. i.e. it should be half the timestep of the simulation /// delta_t = timestep covered by this function. i.e. it should be half the timestep of the simulation
/// feature_size = how many units apart is the center of each adjacent cell on the grid. /// feature_size = how many units apart is the center of each adjacent cell on the grid.
fn step(&mut self, neighbors: Neighbors<'_, R>, delta_t: R, inv_feature_size: R) { fn step(&mut self, neighbors: Neighbors<'_, R>, delta_t: R, inv_feature_size: R) {
@@ -1146,7 +1143,7 @@ pub struct CellState<R> {
h: Vec3<R>, h: Vec3<R>,
} }
impl<R: real::Real> CellState<R> { impl<R: Real> CellState<R> {
pub fn h(&self) -> Vec3<R> { pub fn h(&self) -> Vec3<R> {
self.h self.h
} }
@@ -1281,7 +1278,7 @@ mod test {
do_conv_flt_test(a, b, c, delta_t, &signal, &expected); do_conv_flt_test(a, b, c, delta_t, &signal, &expected);
} }
fn assert_vec3_eq<R: real::Real>(actual: Vec3<R>, expected: Vec3<R>) { fn assert_vec3_eq<R: Real>(actual: Vec3<R>, expected: Vec3<R>) {
let actual = actual.cast::<f64>(); let actual = actual.cast::<f64>();
let expected = expected.cast::<f64>(); let expected = expected.cast::<f64>();
assert_float_eq!(actual.x(), expected.x(), r2nd <= 1e-6); assert_float_eq!(actual.x(), expected.x(), r2nd <= 1e-6);
@@ -1521,7 +1518,7 @@ mod test {
#[test] #[test]
fn accessors() { fn accessors() {
let mut state = StaticSim::new(Index((15, 17, 19).into()), 1e-6); let mut state = SimState::<Static<checked::R64>, checked::R64>::new(Index((15, 17, 19).into()), 1e-6);
assert_eq!(state.width(), 15); assert_eq!(state.width(), 15);
assert_eq!(state.height(), 17); assert_eq!(state.height(), 17);
assert_eq!(state.depth(), 19); assert_eq!(state.depth(), 19);
@@ -1534,7 +1531,7 @@ mod test {
#[test] #[test]
fn energy_conservation_over_time() { fn energy_conservation_over_time() {
let mut state = StaticSim::new(Index((2001, 1, 1).into()), 1e-6); let mut state = SimState::<Static<checked::R64>, checked::R64>::new(Index((2001, 1, 1).into()), 1e-6);
for t in 0..100 { for t in 0..100 {
let angle = (t as f32)/100.0*f32::two_pi(); let angle = (t as f32)/100.0*f32::two_pi();
let sig = angle.sin(); let sig = angle.sin();
@@ -1552,7 +1549,7 @@ mod test {
#[test] #[test]
fn sane_boundary_conditions() { fn sane_boundary_conditions() {
let mut state = StaticSim::new(Index((21, 21, 21).into()), 1e-6); let mut state = SimState::<Static<checked::R64>, checked::R64>::new(Index((21, 21, 21).into()), 1e-6);
for t in 0..10 { for t in 0..10 {
let angle = (t as f32)/10.0*f32::two_pi(); let angle = (t as f32)/10.0*f32::two_pi();
let sig = angle.sin(); let sig = angle.sin();
@@ -1569,8 +1566,8 @@ mod test {
/// Fill the world with the provided material and a stimulus. /// Fill the world with the provided material and a stimulus.
/// Measure energy at the start, and then again after advancing many steps. /// Measure energy at the start, and then again after advancing many steps.
/// Return these two measurements (energy(t=0), energy(t=~=1000)) /// Return these two measurements (energy(t=0), energy(t=~=1000))
fn conductor_test<M: Into<mat::Static<Real>>>(mat: M) -> (f32, f32) { fn conductor_test<M: Into<mat::Static<checked::R64>>>(mat: M) -> (f32, f32) {
let mut state = StaticSim::new(Index((201, 1, 1).into()), 1e-6); let mut state = SimState::<Static<checked::R64>, checked::R64>::new(Index((201, 1, 1).into()), 1e-6);
state.fill_region(&WorldRegion, mat.into()); state.fill_region(&WorldRegion, mat.into());
for t in 0..100 { for t in 0..100 {
let angle = (t as f32)/100.0*f32::two_pi(); let angle = (t as f32)/100.0*f32::two_pi();
@@ -1616,7 +1613,7 @@ mod test {
assert_float_eq!(energy_1/energy_0, 0.0, abs <= 1e-6); assert_float_eq!(energy_1/energy_0, 0.0, abs <= 1e-6);
} }
fn state_for_pml(size: Index) -> SimState<Pml<Real>> { fn state_for_pml(size: Index) -> SimState<Pml<checked::R64>, checked::R64> {
let mut state = SimState::new(size, 1e-6); let mut state = SimState::new(size, 1e-6);
let timestep = state.timestep(); let timestep = state.timestep();
state.fill_boundary_using(size/4, |boundary_ness| { state.fill_boundary_using(size/4, |boundary_ness| {
@@ -1635,15 +1632,13 @@ mod test {
} }
/// Like state_for_pml, but the boundary is an ordinary electric conductor -- not PML /// Like state_for_pml, but the boundary is an ordinary electric conductor -- not PML
fn state_for_graded_conductor(size: Index) -> StaticSim { fn state_for_graded_conductor(size: Index) -> SimState<Conductor<checked::R64>, checked::R64> {
let mut state = StaticSim::new(size, 1e-6); let mut state = SimState::new(size, 1e-6);
let timestep = state.timestep(); let timestep = state.timestep();
state.fill_boundary_using(size/4, |boundary_ness| { state.fill_boundary_using(size/4, |boundary_ness| {
let b = boundary_ness.elem_pow(3.0); let b = boundary_ness.elem_pow(3.0);
let conductivity = Vec3::uniform(f32::eps0() * b.mag() * 0.5 / timestep); let conductivity = Vec3::uniform(f32::eps0() * b.mag() * 0.5 / timestep);
Static { Conductor::new_anisotropic(conductivity)
conductivity: conductivity.cast(), ..Default::default()
}
}); });
state state
} }
@@ -1782,7 +1777,7 @@ mod test {
pml_test_against_baseline(&mut pml_state, &mut baseline_state, Vec3::unit_x()); pml_test_against_baseline(&mut pml_state, &mut baseline_state, Vec3::unit_x());
} }
fn state_for_monodirectional_pml(size: Index) -> SimState<Pml<Real>> { fn state_for_monodirectional_pml(size: Index) -> SimState<Pml<checked::R64>, checked::R64> {
let mut state = SimState::new(size, 1e-6); let mut state = SimState::new(size, 1e-6);
let timestep = state.timestep(); let timestep = state.timestep();
state.fill_boundary_using(size/4, |boundary_ness| { state.fill_boundary_using(size/4, |boundary_ness| {
@@ -1805,7 +1800,7 @@ mod test {
fn pml_ineffective_mono_linear_test<F: Fn(Vec3<f32>) -> Vec3<f32>>( fn pml_ineffective_mono_linear_test<F: Fn(Vec3<f32>) -> Vec3<f32>>(
size: Index, e: Vec3<f32>, shuffle: F size: Index, e: Vec3<f32>, shuffle: F
) { ) {
let mut state = SimState::<Pml<Real>>::new(size, 1e-6); let mut state = SimState::new(size, 1e-6);
let timestep = state.timestep(); let timestep = state.timestep();
state.fill_boundary_using(size/4, |boundary_ness| { state.fill_boundary_using(size/4, |boundary_ness| {
let b = boundary_ness.elem_pow(3.0); let b = boundary_ness.elem_pow(3.0);
@@ -1813,7 +1808,7 @@ mod test {
// let cs = b * 0.5; // let cs = b * 0.5;
// permute the stretching so that it shouldn't interfere with the wave // permute the stretching so that it shouldn't interfere with the wave
let coord_stretch = shuffle(cs); let coord_stretch = shuffle(cs);
Pml::new(coord_stretch) Pml::<checked::R64>::new(coord_stretch)
}); });
let center = Index(state.size().0 / 2); let center = Index(state.size().0 / 2);
let en = pml_test_at(&mut state, e, center); let en = pml_test_at(&mut state, e, center);