diff --git a/crates/coremem/src/bin/bench.rs b/crates/coremem/src/bin/bench.rs index de7eaa1..4852dc1 100644 --- a/crates/coremem/src/bin/bench.rs +++ b/crates/coremem/src/bin/bench.rs @@ -1,5 +1,4 @@ use coremem::{self, Driver, AbstractSim}; -use coremem::sim::legacy::SimState; use coremem::sim::spirv::{SpirvSim, WgpuBackend}; use coremem::sim::units::Frame; use coremem::cross::mat::FullyGenericMaterial; @@ -28,26 +27,14 @@ fn main() { measure_steps("spirv/80", 1, Driver::new( SpirvSim::, WgpuBackend>::new(Index::new(80, 80, 80), 1e-3) )); - measure_steps("sim/80", 1, Driver::<_, SimState>::new( - SimState::new(Index::new(80, 80, 80), 1e-3) - )); measure_steps("spirv/80 step(2)", 2, Driver::new( SpirvSim::, WgpuBackend>::new(Index::new(80, 80, 80), 1e-3) )); - measure_steps("sim/80 step(2)", 2, Driver::<_, SimState>::new( - SimState::new(Index::new(80, 80, 80), 1e-3) - )); measure_steps("spirv/80 step(10)", 10, Driver::new( SpirvSim::, WgpuBackend>::new(Index::new(80, 80, 80), 1e-3) )); - measure_steps("sim/80 step(10)", 10, Driver::<_, SimState>::new( - SimState::new(Index::new(80, 80, 80), 1e-3) - )); measure_steps("spirv/80 step(100)", 100, Driver::new( SpirvSim::, WgpuBackend>::new(Index::new(80, 80, 80), 1e-3) )); - measure_steps("sim/80 step(100)", 100, Driver::<_, SimState>::new( - SimState::new(Index::new(80, 80, 80), 1e-3) - )); } diff --git a/crates/coremem/src/sim/legacy/mat/bh_ferromagnet.rs b/crates/coremem/src/sim/legacy/mat/bh_ferromagnet.rs deleted file mode 100644 index 7809b90..0000000 --- a/crates/coremem/src/sim/legacy/mat/bh_ferromagnet.rs +++ /dev/null @@ -1,365 +0,0 @@ -use super::Material; -use crate::material_compat; -use crate::geom::{Line2d, Polygon2d}; -use crate::real::Real; -use crate::sim::legacy::{CellState, StepParametersMut}; -use crate::cross::vec::{Vec2, Vec3}; - -use lazy_static::lazy_static; -use log::trace; -use serde::{Serialize, Deserialize}; -use std::any::{Any, TypeId}; -use std::cmp::Ordering; -use std::collections::HashMap; -use std::sync::Mutex; - -fn step_linear_ferro(m_mut: &mut Vec3, mh_curve: &MHCurve, context: &CellState, delta_b: Vec3) { - trace!("step_b enter"); - let (h, m) = (context.h(), *m_mut); - let target_hm = h + m + delta_b * R::mu0_inv(); - - // TODO: this is probably not the best way to generalize a BH curve into 3d. - let (_hx, mx) = mh_curve.move_to( - h.x(), - m.x(), - target_hm.x(), - ); - let (_hy, my) = mh_curve.move_to( - h.y(), - m.y(), - target_hm.y(), - ); - let (_hz, mz) = mh_curve.move_to( - h.z(), - m.z(), - target_hm.z(), - ); - - *m_mut = Vec3::new(mx, my, mz); - // let ret = Vec3::new(hx, hy, hz); - trace!("step_b end"); -} - -/// M as a function of H -#[derive(Clone, PartialEq)] -struct MHCurve { - geom: Polygon2d, -} - -#[allow(unused)] -impl MHCurve { - /// Construct a M(H) curve from a sweep from M = 0 to Ms and back down to M = 0. - /// The curve below M = 0 is derived by symmetry. - fn new(points: &[Vec2]) -> Self { - let full_pts: Vec<_> = - points.iter().cloned() - .chain(points.iter().cloned().map(|p| -p)) - .map(|p| p.cast()) - .collect(); - - Self { - geom: Polygon2d::new(full_pts) - } - } - - fn from_bh(points: &[(R2, R2)]) -> Self { - let mh_points: Vec<_> = points.iter().cloned().map(|(h, b)| { - Vec2::new(h, b / R2::mu0() - h) - }).collect(); - - Self::new(&*mh_points) - } - - fn from_mh(points: &[(R2, R2)]) -> Self { - let mh_points: Vec<_> = points.iter().cloned().map(|(h, m)| { - Vec2::new(h, m) - }).collect(); - - Self::new(&*mh_points) - } - - /// Return (Hmax, Mmax) - pub fn extremes(&self) -> Vec2 { - Vec2::new(self.geom.max_x(), self.geom.max_y()) - } - - /// Moves (h, m) towards some location in the MH curve where H + M = target_hm. - /// Returns `Ok((h, m))` if complete; `Err((h, m))` if there's more work to be done (call it - /// again). - fn step_toward(&self, h: R, m: R, target_hm: R) -> Result, Vec2> { - let is_ascending = match target_hm.partial_cmp(&(h + m)).unwrap_or_else(|| panic!("{} {}", h, m)) { - Ordering::Greater => true, - Ordering::Less => false, - _ => return Ok(Vec2::new(h, m)) - }; - if (is_ascending && m == self.geom.max_y()) || (!is_ascending && m == self.geom.min_y()) { - // Fully saturated. m is fixed, while h moves freely - return Ok(Vec2::new(target_hm - m, m)); - } - // Locate the segment which would contain the current point - let mut segments = self.geom.segments(); - let active_segment = loop { - let line = segments.next().unwrap_or_else(|| { - panic!("failed to find segment for h:{}, m:{}, {:?}", h, m, self.geom.segments().collect::>()); - }); - if line.contains_y(m) && line.is_ascending() == is_ascending { - if line.contains_x(h) && line.distance_sq(Vec2::new(h, m)) < R::from_primitive(1.0e-6) { - // (h, m) resides on this line - break line; - } else { - // need to move the point toward this line - let h_intercept = line.x(m); - break Line2d::new(Vec2::new(h, m), Vec2::new(h_intercept, m)); - } - } - }; - trace!("active segment: {:?}", active_segment); - - // Find some m(h) on the active_segment such that sum(h) = h + m(h) = target_hm - let sum_h = active_segment + Line2d::new(Vec2::zero(), Vec2::unit()); - trace!("sum_h: {:?}", sum_h); - let new_h = if sum_h.to().y() != sum_h.from().y() { - sum_h.move_toward_y_unclamped(h, target_hm) - } else { - // avoid a division-by-zero. - // We could be anywhere along this line, but we prefer the endpoint - // so as to escape out of any permanent loops - active_segment.to().x() - }; - trace!("new_h: {}", new_h); - - if sum_h.contains_x(new_h) { - // the segment contains a point with the target H+M - Ok(active_segment.at_x(new_h)) - } else { - // the segment doesn't contain the desired point: clamp and try the next segment - Err(active_segment.clamp_by_x(new_h)) - } - } - fn move_to(&self, mut h: R, mut m: R, target_hm: R) -> (R, R) { - let mut i = 0; - loop { - i += 1; - match self.step_toward(h, m, target_hm) { - Ok(v) => break (v.x(), v.y()), - Err(v) => { - h = v.x(); - m = v.y(); - }, - } - if i % 2048 == 0 { - panic!("unusually high iteration count without converging: {}. args: {}, {}, {}", i, h, m, target_hm); - } - } - } -} - -#[derive(Default, Copy, Clone, PartialEq, Serialize, Deserialize)] -pub struct Ferroxcube3R1 { - m: Vec3, -} - -impl Ferroxcube3R1 { - pub fn new() -> Self { - Self::default() - } -} - -impl Ferroxcube3R1 { - fn curve() -> &'static MHCurve { - lazy_static! { - static ref CURVES: Mutex>> = Mutex::new(HashMap::new()); - } - let mut lock = CURVES.lock().unwrap(); - let curve = lock.entry(TypeId::of::()).or_insert_with(|| { - Box::new(MHCurve::::from_bh(&[ - ( 35.0, 0.0), - ( 50.0, 0.250), - ( 100.0, 0.325), - ( 200.0, 0.350), - (1000.0, 0.390), - // Falling - ( 200.0, 0.360), - ( 100.0, 0.345), - ( 50.0, 0.340), - ( 0.0, 0.325), - ])) - }).downcast_ref::>().unwrap(); - unsafe { std::mem::transmute::<&MHCurve, &'static MHCurve>(curve) } - } -} - -impl Material for Ferroxcube3R1 { - fn step_b(&mut self, context: &CellState, delta_b: Vec3) { - step_linear_ferro(&mut self.m, Self::curve(), context, delta_b) - } - fn m(&self) -> Vec3 { - self.m - } - fn step_parameters_mut<'a>(&'a mut self) -> StepParametersMut<'a, R> { - StepParametersMut::default().with_conductivity(Vec3::uniform(1e-3)) - } -} -material_compat!(R, Ferroxcube3R1); - -/// Simple, square-loop ferrite -#[derive(Default, Copy, Clone, PartialEq, Serialize, Deserialize)] -pub struct MinimalSquare { - m: Vec3, -} - -impl MinimalSquare { - fn curve() -> &'static MHCurve { - lazy_static! { - static ref CURVES: Mutex>> = Mutex::new(HashMap::new()); - } - let mut lock = CURVES.lock().unwrap(); - let curve = lock.entry(TypeId::of::()).or_insert_with(|| { - Box::new(MHCurve::::from_bh(&[ - ( 1.0, 0.0), - ( 2.0, 1000000.0), - // Falling - ( 0.0, 900000.0), - ])) - }).downcast_ref::>().unwrap(); - unsafe { std::mem::transmute::<&MHCurve, &'static MHCurve>(curve) } - } -} - -impl Material for MinimalSquare { - fn step_b(&mut self, context: &CellState, delta_b: Vec3) { - step_linear_ferro(&mut self.m, Self::curve(), context, delta_b) - } - fn m(&self) -> Vec3 { - self.m - } - fn step_parameters_mut<'a>(&'a mut self) -> StepParametersMut<'a, R> { - StepParametersMut::default().with_conductivity(Vec3::uniform(1e-3)) - } -} -material_compat!(R, MinimalSquare); - - -#[cfg(test)] -mod test { - use super::*; - fn mh_curve_for_test() -> MHCurve { - MHCurve::new(&[ - // rising - Vec2::new( 10.0, 0.0), - Vec2::new( 20.0, 100.0), - Vec2::new( 30.0, 150.0), - // falling - Vec2::new( 0.0, 120.0), - // negative rising - Vec2::new(-10.0, 0.0), - Vec2::new(-20.0, -100.0), - Vec2::new(-30.0, -150.0), - // negative falling - Vec2::new( 0.0, -120.0), - ]) - } - - fn assert_step_toward_symmetric(h: f32, m: f32, target_mh: f32, target: Result, Vec2>) { - let curve = mh_curve_for_test(); - let target = match target { - Ok(v) => Ok(v), - Err(v) => Err(v), - }; - let neg_target = match target { - Ok(v) => Ok(-v), - Err(v) => Err(-v), - }; - assert_eq!(curve.step_toward(h, m, target_mh), target); - assert_eq!(curve.step_toward(-h, -m, -target_mh), neg_target); - } - - fn assert_move_to_symmetric(h: f32, m: f32, target_mh: f32, target: (f32, f32)) { - let curve = mh_curve_for_test(); - assert_eq!(curve.move_to(h, m, target_mh), target); - assert_eq!(curve.move_to(-h, -m, -target_mh), (-target.0, -target.1)); - } - - #[test] - fn mh_curve_move_from_inner_to_inner() { - assert_step_toward_symmetric(0.0, 0.0, 5.0, Ok(Vec2::new(5.0, 0.0))); - assert_step_toward_symmetric(0.0, 5.0, 10.0, Ok(Vec2::new(5.0, 5.0))); - - assert_step_toward_symmetric(-5.0, 5.0, -3.0, Ok(Vec2::new(-8.0, 5.0))); - assert_step_toward_symmetric(-5.0, 5.0, 7.0, Ok(Vec2::new(2.0, 5.0))); - - assert_step_toward_symmetric(5.0, -5.0, -3.0, Ok(Vec2::new(2.0, -5.0))); - assert_step_toward_symmetric(5.0, -5.0, 3.0, Ok(Vec2::new(8.0, -5.0))); - } - - #[test] - fn mh_curve_magnetize_along_edge() { - // start of segment NOOP - assert_step_toward_symmetric(10.0, 0.0, 10.0, Ok(Vec2::new(10.0, 0.0))); - // start of segment to middle of segment - assert_step_toward_symmetric(10.0, 0.0, 32.0, Ok(Vec2::new(12.0, 20.0))); - // middle of segment NOOP - assert_step_toward_symmetric(12.0, 20.0, 32.0, Ok(Vec2::new(12.0, 20.0))); - // middle of segment to middle of segment - assert_step_toward_symmetric(12.0, 20.0, 54.0, Ok(Vec2::new(14.0, 40.0))); - // middle of segment to end of segment - assert_step_toward_symmetric(12.0, 20.0, 120.0, Err(Vec2::new(20.0, 100.0))); - } - - #[test] - fn mh_curve_demagnetize_along_edge() { - // start of segment NOOP - assert_step_toward_symmetric(30.0, 150.0, 180.0, Ok(Vec2::new(30.0, 150.0))); - // start of segment to middle of segment - assert_step_toward_symmetric(30.0, 150.0, 160.0, Ok(Vec2::new(20.0, 140.0))); - // middle of segment NOOP - assert_step_toward_symmetric(20.0, 140.0, 160.0, Ok(Vec2::new(20.0, 140.0))); - // middle of segment to middle of segment - assert_step_toward_symmetric(20.0, 140.0, 140.0, Ok(Vec2::new(10.0, 130.0))); - // middle of segment to end of segment - assert_step_toward_symmetric(20.0, 140.0, 120.0, Err(Vec2::new(0.0, 120.0))); - } - - #[test] - fn mh_curve_magnetize_across_edges() { - // Rising from start to middle - assert_move_to_symmetric(10.0, 0.0, 132.0, (22.0, 110.0)); - // Rising from start to saturation - assert_move_to_symmetric(10.0, 0.0, 180.0, (30.0, 150.0)); - // Rising from start to post-saturation - assert_move_to_symmetric(10.0, 0.0, 400.0, (250.0, 150.0)); - // Rising from negative saturation to start - assert_move_to_symmetric(-30.0, -150.0, 10.0, (10.0, 0.0)); - // Rising from negative post-saturation to start - assert_move_to_symmetric(-250.0, -150.0, 10.0, (10.0, 0.0)); - // Rising from negative middle to middle - assert_move_to_symmetric(-22.0, -110.0, 132.0, (22.0, 110.0)); - } - - #[test] - fn mh_curve_demagnetize_across_edges() { - // Falling from saturation to start - assert_move_to_symmetric(30.0, 150.0, 120.0, (0.0, 120.0)); - // Falling from post-saturation to post-saturation - assert_move_to_symmetric(250.0, 150.0, 200.0, (50.0, 150.0)); - // Falling from post-saturation to saturation - assert_move_to_symmetric(250.0, 150.0, 180.0, (30.0, 150.0)); - // Falling from post-saturation to start - assert_move_to_symmetric(250.0, 150.0, 120.0, (0.0, 120.0)); - // Falling from post-saturation to negative saturation - assert_move_to_symmetric(250.0, 150.0, -180.0, (-30.0, -150.0)); - // Falling from post-saturation to negative post-saturation - assert_move_to_symmetric(250.0, 150.0, -400.0, (-250.0, -150.0)); - // Falling from interior to middle - assert_move_to_symmetric(28.0, 130.0, 140.0, (10.0, 130.0)); - // Falling from interior to middle - assert_move_to_symmetric(28.0, 130.0, 130.0, (5.0, 125.0)); - } - - /// Float rounding would cause `inf`s, which manifested as infinite looping. - #[test] - fn regression_no_convergence_3r1() { - let curve = Ferroxcube3R1::curve(); - curve.move_to(-202.04596, -278400.53, -278748.66); - } -} diff --git a/crates/coremem/src/sim/legacy/mat/db.rs b/crates/coremem/src/sim/legacy/mat/db.rs deleted file mode 100644 index 27aaec8..0000000 --- a/crates/coremem/src/sim/legacy/mat/db.rs +++ /dev/null @@ -1,35 +0,0 @@ -//! database of common materials - -use super::{AnisomorphicConductor, IsomorphicConductor, LinearMagnet, Ferroxcube3R1, MinimalSquare}; -use crate::real::Real; -use crate::cross::vec::Vec3; - -pub fn conductor(conductivity: R2) -> IsomorphicConductor { - IsomorphicConductor::new(conductivity.cast()) -} -pub fn anisotropic_conductor(conductivity: Vec3) -> AnisomorphicConductor { - AnisomorphicConductor::new(conductivity) -} - -pub fn copper() -> IsomorphicConductor { - conductor(50_000_000.0) -} - -// See https://en.wikipedia.org/wiki/Permeability_(electromagnetism)#Values_for_some_common_materials -/// This is a simplified form of iron annealed in H. -pub fn linear_annealed_iron() -> LinearMagnet { - LinearMagnet::new(200_000.0) -} -/// This is a simplified form of iron -pub fn linear_iron() -> LinearMagnet { - LinearMagnet::new(5000.0) -} - -/// https://www.ferroxcube.com/upload/media/product/file/MDS/3r1.pdf -pub fn ferroxcube_3r1() -> Ferroxcube3R1 { - Ferroxcube3R1::default() -} -pub fn minimal_square_ferrite() -> MinimalSquare { - MinimalSquare::default() -} - diff --git a/crates/coremem/src/sim/legacy/mat/linear.rs b/crates/coremem/src/sim/legacy/mat/linear.rs deleted file mode 100644 index 5152168..0000000 --- a/crates/coremem/src/sim/legacy/mat/linear.rs +++ /dev/null @@ -1,103 +0,0 @@ -use super::Material; -use crate::material_compat; -use crate::real::Real; -use crate::sim::legacy::CellState; -use crate::cross::vec::Vec3; - -use serde::{Serialize, Deserialize}; - -/// Material which can be magnetized, but has no hysteresis and no coercivity. -#[derive(Copy, Clone, Default, PartialEq, Serialize, Deserialize)] -pub struct LinearMagnet { - /// \mu_r - relative_permeability: Vec3, - m: Vec3, -} - -impl LinearMagnet { - pub fn new(relative_permeability: R2) -> Self { - Self { - relative_permeability: Vec3::uniform(relative_permeability).cast(), - m: Vec3::zero(), - } - } - pub fn new_anisotropic(relative_permeability: Vec3) -> Self { - Self { - relative_permeability: relative_permeability.cast(), - m: Vec3::zero() - } - } -} - -impl Material for LinearMagnet { - fn m(&self) -> Vec3 { - self.m - } - fn step_b(&mut self, _context: &CellState, delta_b: Vec3) { - //```tex - // $B = \mu_0 (H + M) = \mu_0 \mu_r H$ - // $\mu_r H = H + M$ - // $M = (\mu_r - 1) H$ - // $B = \mu_0 (1/(\mu_r - 1) M + M)$ - // $B = \mu_0 \mu_r/(\mu_r - 1) M$ - //``` - let mu_r = self.relative_permeability; - let delta_m = (delta_b*R::mu0_inv()).elem_mul(mu_r - Vec3::unit()).elem_div(mu_r); - self.m += delta_m; - } -} -material_compat!(R, LinearMagnet); - -#[cfg(test)] -mod test { - use super::*; - use float_eq::assert_float_eq; - - #[test] - fn linear_magnet_steep() { - let mut mag = LinearMagnet::::new(5000.0); - - // M = B/mu0 * (mu_r-1)/(mu_r) - mag.step_b(&CellState::default(), Vec3::uniform(1.0)); - assert_float_eq!(mag.m().x(), 795615.56, abs <= 1.0); - - mag.step_b(&CellState::default(), Vec3::uniform(1.0)); - assert_float_eq!(mag.m().x(), 1591231.12, abs <= 1.0); - - mag.step_b(&CellState::default(), Vec3::uniform(-1.0)); - assert_float_eq!(mag.m().x(), 795615.56, abs <= 1.0); - mag.step_b(&CellState::default(), Vec3::uniform(-1.0)); - assert_float_eq!(mag.m().x(), 0.0, abs <= 1.0); - } - - #[test] - fn linear_magnet_shallow() { - let mut mag = LinearMagnet::::new(2.0); - - mag.step_b(&CellState::default(), Vec3::uniform(1.0)); - assert_float_eq!(mag.m().x(), 397887.36, abs <= 1.0); - - mag.step_b(&CellState::default(), Vec3::uniform(-3.0)); - assert_float_eq!(mag.m().x(), -795774.72, abs <= 1.0); - } - - #[test] - fn linear_magnet_accuracy() { - let mut mag = LinearMagnet::::new(5000.0); - - let mut b = Vec3::zero(); - while b.x() < 1.0 { - let delta_b = Vec3::uniform(0.00002); - mag.step_b(&CellState::default(), delta_b); - b += delta_b; - } - while b.x() > 0.0 { - let delta_b = Vec3::uniform(-0.00001); - mag.step_b(&CellState::default(), delta_b); - b += delta_b; - } - // TODO: This error is WAY too big! - // Need to make sure that M+H == mu0*B always - assert_float_eq!(mag.m().x(), b.x() * f32::mu0_inv(), abs <= 900.0); - } -} diff --git a/crates/coremem/src/sim/legacy/mat/mod.rs b/crates/coremem/src/sim/legacy/mat/mod.rs deleted file mode 100644 index d3ffec4..0000000 --- a/crates/coremem/src/sim/legacy/mat/mod.rs +++ /dev/null @@ -1,395 +0,0 @@ -use crate::real::Real; -use crate::sim::legacy::{CellState, PmlParameters, PmlState, StepParameters, StepParametersMut}; -use crate::cross::vec::Vec3; - -use serde::{Serialize, Deserialize}; - -pub mod db; -mod bh_ferromagnet; -mod linear; - -pub use bh_ferromagnet::*; -pub use coremem_cross::mat::{ - AnisomorphicConductor, - Ferroxcube3R1MH, - FullyGenericMaterial, - IsoConductorOr, - IsomorphicConductor, - MHPgram, -}; -pub use linear::*; - -pub trait Material { - fn step_parameters_mut<'a>(&'a mut self) -> StepParametersMut<'a, R> { - // by default, behave as a vacuum - StepParametersMut::default() - } - /// Return the magnetization. - fn m(&self) -> Vec3 { - Vec3::zero() - } - /// Called just before magnetic field is updated. Optionally change any internal state (e.g. magnetization). - fn step_b(&mut self, _context: &CellState, _delta_b: Vec3) { - } -} -#[macro_export] -macro_rules! material_compat { - (R, $mat:path) => { - // XXX this is not that useful an implementation. - // it exists mostly because some users want the `Material::conductivity()` method. - impl crate::mat::Material for $mat { - fn conductivity(&self) -> Vec3 { - crate::sim::legacy::mat::MaterialExt::step_parameters(self).conductivity() - } - fn move_b_vec(&self, _m: Vec3, _target_b: Vec3) -> Vec3 { - unimplemented!() - } - } - } -} - - -pub trait MaterialExt { - fn step_parameters<'a>(&'a self) -> StepParameters<'a, R>; - fn conductivity(&self) -> Vec3; -} - -impl> MaterialExt for M { - fn step_parameters<'a>(&'a self) -> StepParameters<'a, R> { - unsafe { &mut *(self as *const M as *mut M) }.step_parameters_mut().into() - } - fn conductivity(&self) -> Vec3 { - self.step_parameters().conductivity() - } -} - -/// Capable of capturing all field-related information about a material at any -/// snapshot moment-in-time. Useful for serializing state. -#[derive(Clone, Default, PartialEq, Serialize, Deserialize)] -pub struct Static { - pub conductivity: Vec3, - // pub pml: Option<(PmlState, PmlParameters)>, - pub m: Vec3, -} - -impl Static { - pub fn from_material>(m: &M) -> Self { - let p = m.step_parameters(); - Self { - conductivity: p.conductivity(), - // pml: p.pml().map(|(s, p)| (*s, p)), - m: m.m(), - } - } - - // pub fn from_pml(pseudo_conductivity: Vec3) -> Self { - // Self::from_material(&Pml::new(pseudo_conductivity)) - // } -} - -impl Material for Static { - fn step_parameters_mut<'a>(&'a mut self) -> StepParametersMut<'a, R> { - StepParametersMut::new( - self.conductivity, - None, // self.pml.as_mut().map(|(s, p)| (s, *p)), - ) - } - fn m(&self) -> Vec3 { - self.m - } -} -material_compat!(R, Static); - -impl From for Static -where T: Into> -{ - fn from(mat: T) -> Self { - let generic = mat.into(); - Self::from_material(&generic) - } -} - -#[derive(Clone, Default, PartialEq, Serialize, Deserialize)] -pub struct Pml(PmlState, PmlParameters); - -impl Pml { - pub fn new(pseudo_conductivity: Vec3) -> Self { - Self(PmlState::new(), PmlParameters::new(pseudo_conductivity)) - } -} - -impl Material for Pml { - fn step_parameters_mut<'a>(&'a mut self) -> StepParametersMut<'a, R> { - StepParametersMut::default().with_pml(&mut self.0, self.1) - } -} -material_compat!(R, Pml); - - -#[derive(Clone, PartialEq, Serialize, Deserialize)] -pub enum GenericMaterial { - Conductor(AnisomorphicConductor), - LinearMagnet(LinearMagnet), - Pml(Pml), - MBPgram(MBPgram), - Ferroxcube3R1(Ferroxcube3R1), - MinimalSquare(MinimalSquare), -} - -impl Default for GenericMaterial { - fn default() -> Self { - Self::Conductor(Default::default()) - } -} - -impl From> for GenericMaterial { - fn from(inner: AnisomorphicConductor) -> Self { - Self::Conductor(inner) - } -} - -impl From> for GenericMaterial { - fn from(inner: IsomorphicConductor) -> Self { - let iso_r = IsomorphicConductor::new(inner.iso_conductivity().cast::()); - Self::Conductor(iso_r.into()) - } -} - -impl From> for GenericMaterial { - fn from(inner: LinearMagnet) -> Self { - Self::LinearMagnet(inner) - } -} - -impl From> for GenericMaterial { - fn from(inner: Pml) -> Self { - Self::Pml(inner) - } -} - -impl From> for GenericMaterial { - fn from(inner: MBPgram) -> Self { - Self::MBPgram(inner) - } -} - -impl From> for GenericMaterial { - fn from(inner: Ferroxcube3R1) -> Self { - Self::Ferroxcube3R1(inner) - } -} - -impl From> for GenericMaterial { - fn from(inner: MinimalSquare) -> Self { - Self::MinimalSquare(inner) - } -} - -impl Material for GenericMaterial { - fn step_parameters_mut<'a>(&'a mut self) -> StepParametersMut<'a, R> { - use GenericMaterial::*; - match self { - Conductor(inner) => inner.step_parameters_mut(), - LinearMagnet(inner) => inner.step_parameters_mut(), - Pml(inner) => inner.step_parameters_mut(), - MBPgram(inner) => inner.step_parameters_mut(), - Ferroxcube3R1(inner) => inner.step_parameters_mut(), - MinimalSquare(inner) => inner.step_parameters_mut(), - } - } - /// Return the magnetization. - fn m(&self) -> Vec3 { - use GenericMaterial::*; - match self { - Conductor(inner) => inner.m(), - LinearMagnet(inner) => inner.m(), - Pml(inner) => inner.m(), - MBPgram(inner) => inner.m(), - Ferroxcube3R1(inner) => Material::m(inner), - MinimalSquare(inner) => Material::m(inner), - } - } - /// Called just before magnetic field is updated. Optionally change any internal state (e.g. magnetization). - fn step_b(&mut self, context: &CellState, delta_b: Vec3) { - use GenericMaterial::*; - match self { - Conductor(inner) => inner.step_b(context, delta_b), - LinearMagnet(inner) => inner.step_b(context, delta_b), - Pml(inner) => inner.step_b(context, delta_b), - MBPgram(inner) => inner.step_b(context, delta_b), - Ferroxcube3R1(inner) => inner.step_b(context, delta_b), - MinimalSquare(inner) => inner.step_b(context, delta_b), - } - } -} -material_compat!(R, GenericMaterial); - -#[derive(Clone, Serialize, Deserialize)] -pub enum GenericMaterialNoPml { - Conductor(AnisomorphicConductor), - LinearMagnet(LinearMagnet), - MBPgram(MBPgram), - Ferroxcube3R1(Ferroxcube3R1), - MinimalSquare(MinimalSquare), -} - -impl Default for GenericMaterialNoPml { - fn default() -> Self { - AnisomorphicConductor::default().into() - } -} - -impl From> for GenericMaterialNoPml { - fn from(inner: AnisomorphicConductor) -> Self { - Self::Conductor(inner) - } -} - -impl Material for GenericMaterialNoPml { - fn step_parameters_mut<'a>(&'a mut self) -> StepParametersMut<'a, R> { - use GenericMaterialNoPml::*; - match self { - Conductor(inner) => inner.step_parameters_mut(), - LinearMagnet(inner) => inner.step_parameters_mut(), - MBPgram(inner) => inner.step_parameters_mut(), - Ferroxcube3R1(inner) => inner.step_parameters_mut(), - MinimalSquare(inner) => inner.step_parameters_mut(), - } - } - /// Return the magnetization. - fn m(&self) -> Vec3 { - use GenericMaterialNoPml::*; - match self { - Conductor(inner) => inner.m(), - LinearMagnet(inner) => inner.m(), - MBPgram(inner) => inner.m(), - Ferroxcube3R1(inner) => Material::m(inner), - MinimalSquare(inner) => Material::m(inner), - } - } - /// Called just before magnetic field is updated. Optionally change any internal state (e.g. magnetization). - fn step_b(&mut self, context: &CellState, delta_b: Vec3) { - use GenericMaterialNoPml::*; - match self { - Conductor(inner) => inner.step_b(context, delta_b), - LinearMagnet(inner) => inner.step_b(context, delta_b), - MBPgram(inner) => inner.step_b(context, delta_b), - Ferroxcube3R1(inner) => inner.step_b(context, delta_b), - MinimalSquare(inner) => inner.step_b(context, delta_b), - } - } -} -material_compat!(R, GenericMaterialNoPml); - - -/// Materials which have only 1 Vec3. -#[derive(Clone, Serialize, Deserialize)] -pub enum GenericMaterialOneField { - Conductor(AnisomorphicConductor), - Ferroxcube3R1(Ferroxcube3R1), - MinimalSquare(MinimalSquare), -} - -impl Default for GenericMaterialOneField { - fn default() -> Self { - AnisomorphicConductor::default().into() - } -} - -impl From> for GenericMaterialOneField { - fn from(inner: AnisomorphicConductor) -> Self { - Self::Conductor(inner) - } -} - -impl Material for GenericMaterialOneField { - fn step_parameters_mut<'a>(&'a mut self) -> StepParametersMut<'a, R> { - use GenericMaterialOneField::*; - match self { - Conductor(inner) => inner.step_parameters_mut(), - Ferroxcube3R1(inner) => inner.step_parameters_mut(), - MinimalSquare(inner) => inner.step_parameters_mut(), - } - } - /// Return the magnetization. - fn m(&self) -> Vec3 { - use GenericMaterialOneField::*; - match self { - Conductor(inner) => inner.m(), - Ferroxcube3R1(inner) => Material::m(inner), - MinimalSquare(inner) => Material::m(inner), - } - } - /// Called just before magnetic field is updated. Optionally change any internal state (e.g. magnetization). - fn step_b(&mut self, context: &CellState, delta_b: Vec3) { - use GenericMaterialOneField::*; - match self { - Conductor(inner) => inner.step_b(context, delta_b), - Ferroxcube3R1(inner) => inner.step_b(context, delta_b), - MinimalSquare(inner) => inner.step_b(context, delta_b), - } - } -} -material_compat!(R, GenericMaterialOneField); - -// coremem_cross adapters -// TODO: move this to a dedicated file - -/// the coremem_cross Materials are stateless; -/// rather than hold onto their own magnetic fields (for example), the simulation holds that. -/// that's counter to the original cpu-only simulation, in which materials hold their own state. -/// -/// this type adapts any stateless coremem_cross::mat::Material type to be a coremem::mat::Material. -#[derive(Default, Copy, Clone, PartialEq, Serialize, Deserialize)] -pub struct AdaptStateless { - mat: M, - m: Vec3, -} - -impl AdaptStateless { - pub fn into_inner(self) -> M { - self.mat - } -} - -impl From for AdaptStateless { - fn from(mat: M) -> Self { - Self { mat, m: Default::default() } - } -} - -impl> Material for AdaptStateless { - fn m(&self) -> Vec3 { - self.m - } - fn step_parameters_mut<'a>(&'a mut self) -> StepParametersMut<'a, R> { - let c = self.mat.conductivity(); - StepParametersMut::default().with_conductivity(c) - } - fn step_b(&mut self, context: &CellState, delta_b: Vec3) { - let target_b = context.with_m(self.m).b() + delta_b; - self.m = self.mat.move_b_vec(self.m, target_b); - } -} - -// conductors: these are truly stateless -impl Material for AnisomorphicConductor { - fn step_parameters_mut<'a>(&'a mut self) -> StepParametersMut<'a, R> { - let c = coremem_cross::mat::Material::conductivity(self); - StepParametersMut::default().with_conductivity(c) - } -} -impl Material for IsomorphicConductor { - fn step_parameters_mut<'a>(&'a mut self) -> StepParametersMut<'a, R> { - let c = coremem_cross::mat::Material::conductivity(self); - StepParametersMut::default().with_conductivity(c) - } -} - -pub type MBPgram = AdaptStateless>; -impl MBPgram { - pub fn new(b_start: R, b_end: R, m_max: R) -> Self { - coremem_cross::mat::MBPgram::new(b_start, b_end, m_max).into() - } -} -material_compat!(R, MBPgram); diff --git a/crates/coremem/src/sim/legacy/mod.rs b/crates/coremem/src/sim/legacy/mod.rs deleted file mode 100644 index 2d7a810..0000000 --- a/crates/coremem/src/sim/legacy/mod.rs +++ /dev/null @@ -1,1669 +0,0 @@ -pub mod mat; - -use crate::geom::{Coord, Index}; -use crate::cross; -use crate::cross::real::{R32, Real}; -use crate::cross::step::SimMeta; -use crate::cross::vec::{Vec3, Vec3u}; -use crate::sim::{AbstractSim, Fields, GenericSim, StaticSim}; -use crate::stim::Stimulus; - -use mat::{GenericMaterial, Material}; - -use log::trace; -use ndarray::{Array3, Zip}; -use serde::{Serialize, Deserialize}; - -#[derive(Default, Copy, Clone, PartialEq, Serialize, Deserialize)] -pub struct PmlState { - /// Auxiliary fields used when stepping E - aux_e: PmlAux, - /// Auxiliary fields used when stepping H - aux_h: PmlAux, -} - -impl PmlState { - pub fn new() -> Self { - Self::default() - } -} - -#[derive(Default, Copy, Clone, PartialEq, Serialize, Deserialize)] -struct PmlAux { - /// ```tex - /// $Sy^-1 \ast dFz/dy$, where F is either E or H. - /// This is used to advance the X coordinate of the opposite field - /// ``` - sy_conv_dfz_dy: ConvFlt, - sz_conv_dfy_dz: ConvFlt, - /// Used to advance the Y coordinate of the opposite field - sz_conv_dfx_dz: ConvFlt, - sx_conv_dfz_dx: ConvFlt, - /// Used to advance the Z coordinate of the opposite field - sx_conv_dfy_dx: ConvFlt, - sy_conv_dfx_dy: ConvFlt, -} - -#[derive(Default, Copy, Clone, PartialEq, Serialize, Deserialize)] -struct ConvFlt(R); - -impl ConvFlt { - ///```tex - /// 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)$ - /// ``` - #[allow(unused)] - fn step(&mut self, a: R, b: R, c: R, y: R, delta_t: R) -> R { - assert!(c > R::zero(), "ConvFlt::step called without decay: {}", c.to_f64()); - //```tex - // let $v(t) = a \delta(t) \ast y(t) + m(t)$ - // where $m(t) = b exp[-c t] u(t) \ast y(t)$ - // let $k(t) = b exp[-c t] u(t)$, where "k" is for "Kernel". - // $m(t) = \int_0^{t} k(T) y(t - T)dT$ - // $m(t) = \int_0^{\Delta t} k(T) y(t - T)dT + \int_{\Delta t}^t k(T) y(t - T)dT$ - // $m(t) = \int_0^{\Delta t} k(T) y(t - T)dT + \int_{\Delta t}^t b exp[-c T] y(t - T)dT$ - // $m(t) = \int_0^{\Delta t} k(T) y(t - T)dT + \int_0^{t-\Delta t} b exp[-c T] exp[-c \Delta t] y(t - T)dT$ - // $m(t) = \int_0^{\Delta t} k(T) y(t - T)dT + exp[-c \Delta t] v_prev$ - // - // Assume that $y(t)$ is constant for this time-step. - // Looking just at the integral term: - // $b y \int_0^{\Delta t} exp[-c T] dT$ - // $= b y [-1/c](exp[-c \Delta t] - 1)$ - // So - // $m(t) = exp[-c \Delta t] v_prev + -b y/c (exp[-c \Delta t] - 1)$ - //``` - // trace!("pre: {}", self.0); - let e = (-c * delta_t).exp(); - // assert!(self.0 < 1e50 && self.0 > -1e50, "self={}, a={} b={} c={} y={} dt={} e={}", self.0, a, b, c, y, delta_t, e); - self.0 = self.0*e - b*y/c * (e - R::one()); - // trace!("post: {}", self.0); - self.0 + a*y - } - - ///```tex - /// Consider $v(t) = (\delta(t) - \sigma exp[-\sigma t] u(t)) \ast y(t)$ - /// This method will advance from self = $v(t)$ to self = $v(t + \Delta t)$ - /// - /// The user should call this with $e = exp[-\sigma \Delta t]$. - /// This is a more computationally-efficient version of `step` for when $a=1$ and - /// $b = -c$ (i.e. `dc_shift` = 0). - /// ``` - fn step_only_sigma(&mut self, e: R, y: R) -> R { - //```tex - // From the `step` method above, - // $v(t) = a \delta(t) \ast y(t) + m(t)$ - // $m(t) = exp[-c \Delta t] v_prev + -b y/c (exp[-c \Delta t] - 1)$ - // where $a = 1$, $c=\sigma$, $b=-\sigma$. - // So: - // $m(t) = e v_prev + \sigma/\sigma (e - 1) y$ - // and $v(t) = m(t) + y(t)$ - // Simplify this: - // $v(t) = e (v_prev + y)$ - // $m(t) = v(t) - y$ - //``` - let r = e*(self.0 + y); - self.0 = r - y; - r - } -} - - -/// Used to instruct the Sim how to advance the state at a given location. -#[derive(Default, PartialEq)] -pub struct StepParametersMut<'a, R> { - conductivity: Vec3, - pml: Option<(&'a mut PmlState, PmlParameters)>, -} - -impl<'a, R: Real> StepParametersMut<'a, R> { - pub fn new( - conductivity: Vec3, pml: Option<(&'a mut PmlState, PmlParameters)> - ) -> Self { - Self { - conductivity: conductivity.cast(), - pml, - } - } - pub fn with_conductivity(mut self, conductivity: Vec3) -> Self { - self.conductivity = conductivity.cast(); - self - } - pub fn with_pml(mut self, state: &'a mut PmlState, params: PmlParameters) -> Self { - self.pml = Some((state, params)); - self - } -} - -// /// Like `PmlParameters`, but with some additional things parameterized. -// #[derive(Copy, Clone, Default, PartialEq, Serialize, Deserialize)] -// pub struct CpmlParameters { -// /// K; this is similar to \epsilon_r. Try setting this to 1.0 at the interface and -// /// increasing it towards a larger constant at the simulation boundary. -// propagation_constant: Vec3, -// /// sigma; this is the main parameter that should be varied by location within -// /// the sim space -// pseudo_conductivity: Vec3, -// /// a; this combats artifacting from DC waves. Try setting this to a small constant -// /// at the interface and decreasing it to zero at the sim boundary. -// dc_shift: Vec3, -// } -// -// impl CpmlParameters { -// pub fn new(pseudo_conductivity: Vec3, propagation_constant: Vec3, dc_shift: Vec3) -> Self { -// Self { -// pseudo_conductivity, propagation_constant, dc_shift -// } -// } -// } - -/// PML applies a modified formulation of Maxwell's equations which are reflectionless -/// for normally incident waves. This allows the material to behave as a "Perfectly Matched -/// Layer" (PML), dissipating energy AS IF it was a conductor with the provided -/// `pseudo_conductivity`. -#[derive(Copy, Clone, Default, PartialEq, Serialize, Deserialize)] -pub struct PmlParameters { - /// sigma; this is the main parameter that should be varied by location within - /// the sim space. Initialize it to zero at the interface and increase it - /// to a large number at the simulation boundary. - pseudo_conductivity: Vec3, -} - -impl PmlParameters { - pub fn new(pseudo_conductivity: Vec3) -> Self { - Self { - pseudo_conductivity: pseudo_conductivity.cast(), - } - } -} - -/// See `StepParametersMut` -#[derive(PartialEq)] -pub struct StepParameters<'a, R> { - conductivity: Vec3, - pml: Option<(&'a PmlState, PmlParameters)>, -} - -impl<'a, R: Real> StepParameters<'a, R> { - pub fn conductivity(&self) -> Vec3 { - self.conductivity - } - pub fn pml(&self) -> Option<(&'a PmlState, PmlParameters)> { - self.pml - } -} - -impl<'a, R> From> for StepParameters<'a, R> { - fn from(p: StepParametersMut<'a, R>) -> Self { - Self { - conductivity: p.conductivity, - pml: p.pml.map(|(s, p)| (&*s, p)), - } - } -} - -#[derive(Default, Clone, Serialize, Deserialize)] -pub struct SimState> { - pub(super) cells: Array3, - pub(super) e: Array3>, - pub(super) h: Array3>, - feature_size: R, - pub(super) step_no: u64, - timestep: R, -} - -impl SimState { - pub fn size(&self) -> Index { - Index(Vec3u::new( - self.cells.shape()[2] as _, - self.cells.shape()[1] as _, - self.cells.shape()[0] as _, - )) - } - pub fn width(&self) -> u32 { - self.size().x() - } - pub fn height(&self) -> u32 { - self.size().y() - } - pub fn depth(&self) -> u32 { - self.size().z() - } -} - -impl SimState { - pub fn feature_size(&self) -> f32 { - self.feature_size.to_f32() - } - pub fn inv_feature_size(&self) -> f32 { - (R::one() / self.feature_size).to_f32() - } - pub fn timestep(&self) -> f32 { - self.timestep.to_f32() - } -} - -impl SimState { - pub fn new(size: Index, feature_size: f32) -> Self { - // 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 - // of a grid. The distance traveled during that time is \sqrt{3}. Hence each timestep needs - // to be \sqrt{3}/3, or 0.577. - // In 2d, this would be \sqrt{2}/2. - // It's an upper limit though; it should be safe to go lower. - let courant = 0.577; - let feature_size = R::from_primitive(feature_size); - let timestep = feature_size / R::c() * R::from_primitive(courant); - Self { - cells: Array3::default((size.z() as _, size.y() as _, size.x() as _)), - e: Array3::default((size.z() as _, size.y() as _, size.x() as _)), - h: Array3::default((size.z() as _, size.y() as _, size.x() as _)), - feature_size, - timestep, - ..Default::default() - } - } -} - -impl + cross::mat::Material + Send + Sync> SimState { - pub fn step(&mut self) { - let time_step = self.timestep; - let half_time_step = R::half() * time_step; - let feature_size = self.feature_size; - let inv_feature_size = R::one() / feature_size; - - // advance all the electric fields - trace!("step_e begin"); - let h_field = &self.h; - Zip::from(ndarray::indices_of(&self.cells)).and(&mut self.cells).and(&mut self.e).par_for_each( - |(z, y, x), mat, e| { - let neighbors = Neighbors::new(h_field, [z, y, x]); - let h = &h_field[[z, y, x]]; - ECell { mat, e, h }.step(neighbors, half_time_step, inv_feature_size); - }); - trace!("step_e end"); - - // advance all the magnetic fields - trace!("step_h begin"); - let e_field = &self.e; - Zip::from(ndarray::indices_of(&self.cells)).and(&mut self.cells).and(&mut self.h).par_for_each( - |(z, y, x), mat, h| { - let neighbors = Neighbors::new(e_field, [z, y, x]); - let e = &e_field[[z, y, x]]; - HCell { mat, h, e }.step(neighbors, half_time_step, inv_feature_size); - }); - trace!("step_h end"); - - self.step_no += 1; - } - - pub fn apply_stimulus + ?Sized>(&mut self, stim: &S) { - let feature_size = self.feature_size().cast::(); - - let t_sec = self.time().cast::(); - let timestep = self.timestep.cast::(); - 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 idx = Index::new(x as u32, y as u32, z as u32); - let densities = stim.at(t_sec, feature_size, idx) * timestep; - *e += densities.e(); - *h += densities.h(); - }); - trace!("apply_stimulus end"); - } -} - -impl + cross::mat::Material + Send + Sync> AbstractSim for SimState { - type Real = R; - type Material = M; - fn meta(&self) -> SimMeta { - let dim = Vec3u::new( - self.cells.shape()[2] as _, - self.cells.shape()[1] as _, - self.cells.shape()[0] as _, - ); - let feature_size = self.feature_size.to_f32(); - let inv_feature_size = 1.0/feature_size; - let time_step = self.timestep.to_f32(); - SimMeta { - dim, feature_size, inv_feature_size, time_step - } - } - fn step_no(&self) -> u64 { - self.step_no - } - - fn put_material_index(&mut self, pos: Index, mat: M) { - *self.get_mat_mut(pos) = mat; - } - - fn get_material_index(&self, pos: Index) -> &M { - self.get_material(pos) - } - - fn step_multiple>(&mut self, num_steps: u32, s: &S) { - for _ in 0..num_steps { - self.apply_stimulus(s); - SimState::step(self); - } - } - - fn fields_at_index(&self, pos: Index) -> Fields { - let idx = [pos.z() as usize, pos.y() as _, pos.x() as _]; - match self.cells.get(idx) { - Some(mat) => Fields::new( - self.e[idx], self.h[idx], mat.m() - ).cast(), - None => Default::default(), - } - } - - fn to_static(&self) -> StaticSim { - unimplemented!() - } - fn to_generic(&self) -> GenericSim { - unimplemented!() - } -} - -#[allow(unused)] -impl SimState { - fn get_material(&self, c: C) -> &M { - let at = c.to_index(self.feature_size.cast()); - &self.cells[at.row_major_idx()] - } - fn get_mat_mut(&mut self, c: C) -> &mut M { - let at = c.to_index(self.feature_size.cast()); - &mut self.cells[at.row_major_idx()] - } - - fn get_e(&self, c: C) -> &Vec3 { - let at = c.to_index(self.feature_size.cast()); - &self.e[at.row_major_idx()] - } - fn get_e_mut(&mut self, c: C) -> &mut Vec3 { - let at = c.to_index(self.feature_size.cast()); - &mut self.e[at.row_major_idx()] - } - - fn get_h(&self, c: C) -> &Vec3 { - let at = c.to_index(self.feature_size.cast()); - &self.h[at.row_major_idx()] - } - fn get_h_mut(&mut self, c: C) -> &mut Vec3 { - let at = c.to_index(self.feature_size.cast()); - &mut self.h[at.row_major_idx()] - } - - fn get_hcell<'a, C: Coord>(&'a mut self, c: C) -> HCell<'a, R, M> { - let at = c.to_index(self.feature_size.cast()); - let idx = at.row_major_idx(); - HCell { - e: &self.e[idx], - h: &mut self.h[idx], - mat: &mut self.cells[idx], - } - } -} - - -struct Neighbors<'a, R> { - left: Option<&'a Vec3>, - right: Option<&'a Vec3>, - up: Option<&'a Vec3>, - down: Option<&'a Vec3>, - out: Option<&'a Vec3>, - in_: Option<&'a Vec3>, -} - -impl<'a, R> Neighbors<'a, R> { - fn new(array: &'a Array3>, location: [usize; 3]) -> Self { - let [z, y, x] = location; - let left = array.get([z, y, x.wrapping_sub(1)]); - let right = array.get([z, y, x + 1]); - let up = array.get([z, y.wrapping_sub(1), x]); - let down = array.get([z, y + 1, x]); - let out = array.get([z.wrapping_sub(1), y, x]); - let in_ = array.get([z + 1, y, x]); - Self { - left, right, up, down, out, 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 -#[allow(unused)] -fn solve_step_diff_eq( - nabla_g: Vec3, - df_dt_coeff: Vec3, - f_coeff: Vec3, - f_int_coeff: Vec3, - f_int_int_coeff: Vec3, - delta_t: R, - f_prev: Vec3, - f_int_prev: Vec3, - f_int_int_prev: Vec3, -) -> Vec3 { - // ```tex - // $f(t)$ from $(T-\Delta t)$ to $(T+\Delta t)$ is a line. As such: - // $f(t') = f_prev + t'*db/dt$ where $t'$ is the time since the last sample of $f$, i.e. $f_prev$ - // and - // $\int f|t' = f_int_prev + t'*f_prev + 1/2 t'^2 df/dt$ - // $\int \int f|t' = f_int_int_prev + t'*f_int_prev + 1/2 t'^2*f_prev + 1/6 t'^3*df/dt$ - // ``` - // Solve coeffs in: - // nabla_g = df_dt_coeff * df/dt + f_prev_coeff * f_prev - // + f_int_prev_coeff * f_int_prev + f_int_int_prev_coeff * f_int_int_prev - let f_int_int_prev_coeff = f_int_int_coeff; - let f_int_prev_coeff = f_int_coeff + f_int_int_coeff*delta_t; - let f_prev_coeff = f_coeff + (f_int_coeff + - f_int_int_coeff*R::half() - )*delta_t; - let df_dt_coeff = df_dt_coeff + (f_coeff + - (f_int_coeff + f_int_int_coeff*R::third()*delta_t)*R::half()*delta_t - )*delta_t; - - // Solve coeffs in nabla_g = df_dt_coeff + offset; - let offset = f_prev_coeff.elem_mul(f_prev) + - f_int_prev_coeff.elem_mul(f_int_prev) + - f_int_int_prev_coeff.elem_mul(f_int_int_prev); - (nabla_g - offset).elem_div(df_dt_coeff) -} - -/// Solve df/dt in: -/// nabla_g = df_dt_coeff * df/dt + f_coeff * f -/// Lower-order (more efficient) version of `solve_step_diff_eq`. -fn solve_step_diff_eq_linear( - nabla_g: Vec3, - df_dt_coeff: Vec3, - f_coeff: Vec3, - delta_t: R, - f_prev: Vec3, -) -> Vec3 { - // Solve coeffs in: - // nabla_g = df_dt_coeff * df/dt + f_prev_coeff * f_prev - let df_dt_coeff = df_dt_coeff + f_coeff*delta_t; - - // Solve coeffs in nabla_g = df_dt_coeff + offset; - let offset = f_coeff.elem_mul(f_prev); - (nabla_g - offset).elem_div(df_dt_coeff) -} - -struct HCell<'a, R, M> { - mat: &'a mut M, - h: &'a mut Vec3, - e: &'a Vec3, -} - -struct ECell<'a, R, M> { - mat: &'a mut M, - e: &'a mut Vec3, - h: &'a Vec3, -} - -impl<'a, R: Real, M: Material> HCell<'a, R, M> { - pub fn b(&self) -> Vec3 { - (*self.h + self.mat.m()) * R::mu0() - } - - fn step(&mut self, neighbors: Neighbors<'_, R>, delta_t: R, inv_feature_size: R) { - // ```tex - // Maxwell-Faraday equation: $\nabla x E = -dB/dt$ (1) - // - // Note: $\nabla x E$ evaluates to: - // | i j k | [ dEz/dy - dEy/dz ] - // | | | | - // | d/dx d/dy d/dz | = | dEx/dz - dEz/dx | - // | | | | - // | Ex Ey Ez | [ dEy/dx - dEx/dy ] - // - // To accomodate PML materials, we allow coordinate stretching, specifically - // in the complex domain. - // - // For the rest of this analysis, we'll use $\omega$. This comes from a wave of the form - // $E(t) = A exp(-j\omega t)$ (2) - // Maxwell's equations are all linear, so we'll do the derivations considering just a - // single such wave, knowing that what we derive fully generalizes to any $E(t)$. - // If we do it right, then $\omega$ will drop out from the end result. - // - // Consider: $x' = x + jf(x)$. - // Then $dx' = (1 + jdf/dx)dx$ - // Let $df/dx = \gamma_x/\omega$. This is chosen to eliminate $\omega$ in the end result. - // Then, $d/dx = 1/(1 + j\gamma_x/\omega) d/dx'$ - // Likewise, $d/dy = 1/(1 + j\gamma_y/\omega) d/dy'$ - // $d/dz = 1/(1 + j\gamma_z/\omega) d/dz'$ - // - // Let $S = 1 + j\gamma/\omega$, such that $d/dx = Sx^-1 d/dx'$, etc. - // - // Then (1) can be expanded to the following matrix multiplication - // $ [ 0 -Sz^-1 d/dz' Sy^-1 d/dy' ] $ - // $ | Sz^-1 d/dz' 0 -Sx^-1 d/dx' | E = -dB/dt$ - // $ [ -Sy^-1 d/dy' Sx^-1 d/dx' 0 ]$ - // - // Normalize by multiplying each side by a diagonal matrix: - // - // $ [ 0 -Sy d/dz' Sz d/dy' ] [ SySz ] $ - // $ -| Sx d/dz' 0 -Sz d/dx' | E = diag| SxSz | dB/dt (3)$ - // $ [ -Sx d/dy' Sy d/dx' 0 ] [ SxSy ] $ - // - // The LHS simplifies to: - // $ -\nabla x (diag(S) E)$ - // $ = -\nabla x (E + diag(\gamma) j/\omega E)$ - // - // From (2), we know $B$ is of the form $B = A exp(-j\omega t)$ - // Therefore $\int_0^t B(\tau) d\tau = j/\omega [B(t) - B(0)] = j/\omega B(t)$ (we assume $B(0) = 0$) - // Thus the LHS becomes: - // $ = -\nabla x (E + diag(\gamma) \int E)$ - // - // Expand (3), letting $I = j/\omega$: - // $ |(1 + \gamma_yI)(1 + \gamma_zI)|$ - // $ -\nabla x (E + diag(\gamma) \int E) = diag[(1 + \gamma_xI)(1 + \gamma_zI)] dB/dt$ - // $ |(1 + \gamma_xI)(1 + \gamma_yI)|$ - // - // $ |\gamma_y + \gamma_z| |\gamma_y\gamma_z|$ - // $ -\nabla x (E + diag(\gamma) \int E) = dB/dt + diag[\gamma_x + \gamma_z] I dB/dt + diag[\gamma_x\gamma_z] I^2 dB/dt$ - // $ |\gamma_x + \gamma_y| |\gamma_x\gamma_y|$ - // - // As seen above, $I = j/w$ is equivalent to integration: - // $ |\gamma_y + \gamma_z| |\gamma_y\gamma_z|$ - // $ -\nabla x (E + diag(\gamma) \int E) = dB/dt + diag[\gamma_x + \gamma_z] B + diag[\gamma_x\gamma_z] \int B$ - // $ |\gamma_x + \gamma_y| |\gamma_x\gamma_y|$ - // - // Finally, discretization. We're expanding about $t = T$, where we know $B(T-\Delta t)$ and $E(T)$. - // $B(t)$ from $(T-\Delta t)$ to $(T+\Delta t)$ is a line. As such: - // $B(t') = B_prev + t'*db/dt$ where $t'$ is the time since the last sample of $B$, i.e. $B_prev$ - // and - // $\int B|t' = B_int_prev + t'*B_prev + 1/2 t'^2 dB/dt - // - // ``` - // let inv_feature_size = Real::from_inner(1.0) / feature_size; - let twice_delta_t = R::two() * delta_t; - - let b_prev = self.b(); - - let (delta_ey_x, delta_ez_x) = match (neighbors.right, neighbors.left) { - (Some(right), _) => ( - inv_feature_size * (right.y() - self.e.y()), - inv_feature_size * (right.z() - self.e.z()), - ), - _ => (R::zero(), R::zero()), - }; - let (delta_ex_y, delta_ez_y) = match (neighbors.down, neighbors.up) { - (Some(down), _) => ( - inv_feature_size * (down.x() - self.e.x()), - inv_feature_size * (down.z() - self.e.z()), - ), - _ => (R::zero(), R::zero()), - }; - let (delta_ex_z, delta_ey_z) = match (neighbors.in_, neighbors.out) { - (Some(in_), _) => ( - inv_feature_size * (in_.x() - self.e.x()), - inv_feature_size * (in_.y() - self.e.y()), - ), - _ => (R::zero(), R::zero()), - }; - - let StepParametersMut { conductivity: _, pml } = self.mat.step_parameters_mut(); - // XXX conductivity is unused when stepping h. - let nabla_e = match pml { - None => { - // \nabla x E - Vec3::from((delta_ez_y - delta_ey_z, delta_ex_z - delta_ez_x, delta_ey_x - delta_ex_y)) - } - Some((pml, PmlParameters { pseudo_conductivity })) => { - //```tex - // S = propagation_const + pseudo_conductivity / (dc_shift + j\omega \mu_0) - // $S = 1 + sigma/(j\omega \mu_0)$ - // After some derivation, we get: - // $F{S^-1} = \delta(t) - \sigma/\mu_0exp[-t(\sigma/\mu_0)]u(t)$ - //``` - let sigma = pseudo_conductivity; - - let e = (-sigma * twice_delta_t).exp(); - let aux = &mut pml.aux_h; - // (S^-1 \nabla) x E - Vec3::from(( - aux.sy_conv_dfz_dy.step_only_sigma(e.y(), delta_ez_y) - - aux.sz_conv_dfy_dz.step_only_sigma(e.z(), delta_ey_z), - aux.sz_conv_dfx_dz.step_only_sigma(e.z(), delta_ex_z) - - aux.sx_conv_dfz_dx.step_only_sigma(e.x(), delta_ez_x), - aux.sx_conv_dfy_dx.step_only_sigma(e.x(), delta_ey_x) - - aux.sy_conv_dfx_dy.step_only_sigma(e.y(), delta_ex_y), - )) - } - }; - - // From: $\nabla x E = -dB/dt$ - let db_dt = -nabla_e; - // let db_dt = solve_step_diff_eq( - // nabla_e, - // -Vec3::unit(), // dB/dt coeff - // Vec3::zero(), // B coeff - // Vec3::zero(), // \int B coeff - // Vec3::zero(), // \int \int B coeff - // delta_t, - // b_prev, - // Vec3::zero(), // b_int_prev, - // Vec3::zero(), // b_int_int_prev - // ); - - // Maxwell's equation gave us delta_b. Material parameters let us turn that into h. - let delta_b = db_dt * twice_delta_t; - let state = CellState { - e: *self.e, - h: *self.h, - }; - self.mat.step_b(&state, delta_b); - // mu0 (H + M) = B - // therefore, H = B mu0^-1 - M - // println!("cpu-step_h: nabla_e: {}", nabla_e); - // println!("cpu-step_h: delta_h: {}", delta_b * R::mu0_inv()); - *self.h = (b_prev + delta_b) * R::mu0_inv() - self.mat.m(); - } -} - -impl<'a, R: Real, M: Material> ECell<'a, R, M> { - /// 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. - fn step(&mut self, neighbors: Neighbors<'_, R>, delta_t: R, inv_feature_size: R) { - // ```tex - // Ampere's circuital law with Maxwell's addition, in SI units ("macroscopic version"): - // $\nabla x H = J_f + dD/dt$ where $J_f$ = current density = $\sigma E$, $\sigma$ being a material parameter ("conductivity") - // Note that $D = \epsilon_0 E + P$, but we don't simulate any material where $P \ne 0$. - // Expand $D$: $\nabla x H = J_f + \epsilon_0 dE/dt$ - // Expand $J$: $\nabla x H = \sigma E + \epsilon_0 dE/dt$ - // Rearrange: $dE/dt = 1/\epsilon_0 (\nabla x H - \sigma E)$ - // - // let $E_p$ be $E$ at $T-\Delta t$ and $E_n$ be $E$ at $T+\Delta t$. - // $(E_n - E_p)/(2\Delta{t}) = 1/\epsilon_0 (\nabla x H - \sigma (E_n + E_p)/2)$ - // Normalize: $E_n - E_p = 2\Delta{t}/\epsilon_0 (\nabla x H - \sigma (E_n + E_p)/2)$ - // Expand: $E_n - E_p = 2\Delta{t}/\epsilon_0 \nabla x H - \sigma \Delta{t}/\epsilon_0 (E_n + E_p)$ - // Rearrange: $E_n(1 + \sigma \Delta{t}/\epsilon_0) = E_p(1 - \sigma \Delta{t}/\epsilon_0) + 2\Delta{t}/\epsilon_0 \nabla x H$ - // Then $E_n$ (i.e. the E vector after this step) is trivially solved - // - // Note: $\nabla x H$ evaluates to: - // | i j k | [ dHz/dy - dHy/dz ] - // | | | | - // | d/dx d/dy d/dz | = | dHx/dz - dHz/dx | - // | | | | - // | Hx Hy Hz | [ dHy/dx - dHx/dy ] - // - // Boundary conditions (TODO): - // For something on the yz plane, E(x=0, t+1) should be E(x=1, t). i.e. the wave travels - // $\Delta x$ = 1 per timestep. Roughly, we want to take the line from E(x=0, t), E(x=1, t) - // and continue that to E(x=-1, t) ? - // - // ``` - - let twice_delta_t = R::two() * delta_t; - - let e_prev = *self.e; - - let (delta_hy_x, delta_hz_x) = match (neighbors.left, neighbors.right) { - (Some(left), _) => ( - inv_feature_size * (self.h.y() - left.y()), - inv_feature_size * (self.h.z() - left.z()), - ), - _ => (R::zero(), R::zero()), - }; - let (delta_hx_y, delta_hz_y) = match (neighbors.up, neighbors.down) { - (Some(up), _) => ( - inv_feature_size * (self.h.x() - up.x()), - inv_feature_size * (self.h.z() - up.z()), - ), - _ => (R::zero(), R::zero()), - }; - let (delta_hx_z, delta_hy_z) = match (neighbors.out, neighbors.in_) { - (Some(out), _) => ( - inv_feature_size * (self.h.x() - out.x()), - inv_feature_size * (self.h.y() - out.y()), - ), - _ => (R::zero(), R::zero()), - }; - - let StepParametersMut { conductivity, pml } = self.mat.step_parameters_mut(); - let nabla_h = match pml { - None => { - // \nabla x H - Vec3::from((delta_hz_y - delta_hy_z, delta_hx_z - delta_hz_x, delta_hy_x - delta_hx_y)) - } - Some((pml, PmlParameters { pseudo_conductivity })) => { - //```tex - // S = propagation_const + pseudo_conductivity / (dc_shift + j\omega \epsilon_0) - // $S = 1 + sigma/(j\omega \epsilon_0)$ - // After some derivation, we get: - // $F{S^-1} = \delta(t) - \sigma/\epsilon_0exp[-t(\sigma/\epsilon_0)]u(t)$ - //``` - let sigma = pseudo_conductivity; - - let e = (-sigma * twice_delta_t).exp(); - let aux = &mut pml.aux_e; - // (S^-1 \nabla) x H - Vec3::from(( - aux.sy_conv_dfz_dy.step_only_sigma(e.y(), delta_hz_y) - - aux.sz_conv_dfy_dz.step_only_sigma(e.z(), delta_hy_z), - aux.sz_conv_dfx_dz.step_only_sigma(e.z(), delta_hx_z) - - aux.sx_conv_dfz_dx.step_only_sigma(e.x(), delta_hz_x), - aux.sx_conv_dfy_dx.step_only_sigma(e.x(), delta_hy_x) - - aux.sy_conv_dfx_dy.step_only_sigma(e.y(), delta_hx_y), - )) - } - }; - - // From: $\nabla x H = \epsilon_0 dE/dt + \sigma E$ - let de_dt = solve_step_diff_eq_linear( - nabla_h, - Vec3::uniform(R::eps0()), // dE/dt coeff - conductivity, // E coeff - delta_t, - e_prev, - ); - - // println!("cpu-step_e: nabla_h: {}", nabla_h); - // println!("cpu-step_e: delta_e: {}", de_dt*twice_delta_t); - *self.e = e_prev + de_dt*twice_delta_t; - } -} - -#[derive(Copy, Clone, Debug, Default, PartialEq, Serialize, Deserialize)] -pub struct CellState { - e: Vec3, - h: Vec3, -} - -impl CellState { - pub fn e(&self) -> Vec3 { - self.e - } - pub fn h(&self) -> Vec3 { - self.h - } - pub fn with_m(&self, m: Vec3) -> Fields { - Fields::new(self.e, self.h, m) - } -} - - -/// Most of these serve as "smoke tests", to catch and pinpoint the most egregious of errors, or to -/// understand the impact of future changes (like changing default boundary conditions). -#[cfg(test)] -mod test { - use super::*; - 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 _}; - use float_eq::assert_float_eq; - use more_asserts::*; - use rand::rngs::StdRng; - use rand::{Rng as _, SeedableRng as _}; - use std::fmt::Debug; - use std::ops::RangeBounds; - - fn do_conv_flt_test(a: f64, b: f64, c: f64, delta_t: f64, signal: &[f64], expected: &[f64]) { - let mut cf = ConvFlt(0.0); - for (s, e) in signal.iter().zip(expected.iter()) { - let v = cf.step( - a, - b, - c, - *s, - delta_t, - ); - assert_float_eq!(v, e, r2nd <= 1e-6); - } - } - - // PORT-STATUS: skip - #[test] - fn conv_flt_trivial() { - let signal = [2.0, 7.0, -3.0]; - // kernel: 1.0 \delta(t) - let expected = [2.0, 7.0, -3.0]; - let (a, b, c, delta_t) = (1.0, 0.0, 1.0, 0.1); - do_conv_flt_test(a, b, c, delta_t, &signal, &expected); - } - - // PORT-STATUS: skip - #[test] - fn conv_flt_trivial_scaled() { - let signal = [2.0, 7.0, -3.0]; - // kernel: 10.0 \delta(t) - let expected = [20.0, 70.0, -30.0]; - let (a, b, c, delta_t) = (10.0, 0.0, 1.0, 0.1); - do_conv_flt_test(a, b, c, delta_t, &signal, &expected); - } - - // PORT-STATUS: skip - #[test] - fn conv_flt_exp_trivial() { - let signal = [2.0, 0.0, 0.0]; - // kernel: e(-t) - // \int_0^1 e(-t) dt = [1 - exp(-1)] - let exp_neg_0 = (-0.0f64).exp(); - let exp_neg_1 = (-1.0f64).exp(); - let exp_neg_2 = (-2.0f64).exp(); - let exp_neg_3 = (-3.0f64).exp(); - let expected = [ - 2.0*(exp_neg_0 - exp_neg_1), - 2.0*(exp_neg_1 - exp_neg_2), - 2.0*(exp_neg_2 - exp_neg_3), - ]; - let (a, b, c, delta_t) = (0.0, 1.0, 1.0, 1.0); - do_conv_flt_test(a, b, c, delta_t, &signal, &expected); - } - - // PORT-STATUS: skip - #[test] - fn conv_flt_exp_time_scaled() { - let signal = [2.0, 0.0, 0.0]; - // kernel: e(-3*t) - // \int_0^0.2 e(-3*t) dt = [1 - exp(-0.6)]/3 - let exp_neg_00 = (-0.0f64).exp(); - let exp_neg_06 = (-0.6f64).exp(); - let exp_neg_12 = (-1.2f64).exp(); - let exp_neg_18 = (-1.8f64).exp(); - let expected = [ - 2.0/3.0*(exp_neg_00 - exp_neg_06), - 2.0/3.0*(exp_neg_06 - exp_neg_12), - 2.0/3.0*(exp_neg_12 - exp_neg_18), - ]; - let (a, b, c, delta_t) = (0.0, 1.0, 3.0, 0.2); - do_conv_flt_test(a, b, c, delta_t, &signal, &expected); - } - - // PORT-STATUS: skip - #[test] - fn conv_flt_exp_multi_input() { - let signal = [2.0, 7.0, -3.0]; - // kernel: e(-3*t) - // \int_0^0.2 e(-3*t) dt = [1 - exp(-0.6)]/3 - let exp_neg_00 = (-0.0f64).exp(); - let exp_neg_06 = (-0.6f64).exp(); - let exp_neg_12 = (-1.2f64).exp(); - let exp_neg_18 = (-1.8f64).exp(); - let expected = [ - 2.0/3.0*(exp_neg_00 - exp_neg_06), - 7.0/3.0*(exp_neg_00 - exp_neg_06) + 2.0/3.0*(exp_neg_06 - exp_neg_12), - -3.0/3.0*(exp_neg_00 - exp_neg_06) - + 7.0/3.0*(exp_neg_06 - exp_neg_12) - + 2.0/3.0*(exp_neg_12 - exp_neg_18), - ]; - let (a, b, c, delta_t) = (0.0, 1.0, 3.0, 0.2); - do_conv_flt_test(a, b, c, delta_t, &signal, &expected); - } - - // PORT-STATUS: skip - #[test] - fn conv_flt_all_params() { - let signal = [2.0, 7.0, -3.0]; - // kernel: 0.3 \delta(t) + 0.4 * e(-3*t) - // \int_0^0.2 e(-3*t) dt = [1 - exp(-0.6)]/3 - let exp_neg_00 = (-0.0f64).exp(); - let exp_neg_06 = (-0.6f64).exp(); - let exp_neg_12 = (-1.2f64).exp(); - let exp_neg_18 = (-1.8f64).exp(); - let expected_exp = [ - 2.0/3.0*(exp_neg_00 - exp_neg_06), - 7.0/3.0*(exp_neg_00 - exp_neg_06) + 2.0/3.0*(exp_neg_06 - exp_neg_12), - -3.0/3.0*(exp_neg_00 - exp_neg_06) - + 7.0/3.0*(exp_neg_06 - exp_neg_12) - + 2.0/3.0*(exp_neg_12 - exp_neg_18), - ]; - let expected = [ - 0.3 * signal[0] + 0.4 * expected_exp[0], - 0.3 * signal[1] + 0.4 * expected_exp[1], - 0.3 * signal[2] + 0.4 * expected_exp[2], - ]; - let (a, b, c, delta_t) = (0.3, 0.4, 3.0, 0.2); - do_conv_flt_test(a, b, c, delta_t, &signal, &expected); - } - - fn assert_vec3_eq(actual: Vec3, expected: Vec3) { - let actual = actual.cast::(); - let expected = expected.cast::(); - assert_float_eq!(actual.x(), expected.x(), r2nd <= 1e-6); - assert_float_eq!(actual.y(), expected.y(), r2nd <= 1e-6); - assert_float_eq!(actual.z(), expected.z(), r2nd <= 1e-6); - } - - // PORT-STATUS: skip - #[test] - fn step_diff_eq_trivial() { - // f(t) = 0.0 - // nabla_g = 0.1 df/dt - let df_dt = solve_step_diff_eq( - Vec3::new(2.0, 3.0, 4.0), // nabla_g - Vec3::uniform(0.1), // df_dt_coeff - Vec3::zero(), - Vec3::zero(), - Vec3::zero(), - 0.5, // delta_t - Vec3::zero(), - Vec3::zero(), - Vec3::zero(), - ); - assert_vec3_eq(df_dt, Vec3::new(20.0, 30.0, 40.0)); - } - - // PORT-STATUS: skip - #[test] - fn step_diff_eq_trivial_3rd_order() { - // f(t) = 0.0 - // nabla_g = 0.1 df/dt - f - 3 f_int + 0.4 f_int_int - let df_dt = solve_step_diff_eq( - Vec3::new(2.0, 3.0, 4.0), // nabla_g - Vec3::uniform(0.1), // df_dt_coeff - Vec3::uniform(-1.0), // f_coeff - Vec3::uniform(-3.0), // f_int_coeff - Vec3::uniform(0.4), // f_int_int_coeff - 0.0, // delta_t - Vec3::zero(), - Vec3::zero(), - Vec3::zero(), - ); - assert_vec3_eq(df_dt, Vec3::new(20.0, 30.0, 40.0)); - } - - // PORT-STATUS: skip - #[test] - fn step_diff_eq_1st_order_no_dt() { - // nabla_g = 0.1 df/dt + 2 f - // Let f(t) = sin(2 t) - // then f'(t) = 2 cos(2 t) - let f = Vec3::uniform((5.0f64).sin()); - let df_dt = Vec3::uniform(2.0 * (5.0f64).cos()); - let nabla_g = df_dt*0.1 + f*2.0; - let actual_df_dt = solve_step_diff_eq( - nabla_g, - Vec3::uniform(0.1), - Vec3::uniform(2.0), - Vec3::zero(), - Vec3::zero(), - 0.0, // delta_t - f, // f_prev - Vec3::zero(), - Vec3::zero(), - ); - assert_vec3_eq(actual_df_dt, df_dt); - } - - // PORT-STATUS: skip - #[test] - fn step_diff_eq_2nd_order_no_dt() { - // nabla_g = 0.1 df/dt + 2 f + 0.4 f_int - // Let f(t) = sin(2 t) - // then f'(t) = 2 cos(2 t) - // then \int f(t) = -1/2 cos(2 t) - let f = Vec3::uniform((5.0f64).sin()); - let df_dt = Vec3::uniform(2.0 * (5.0f64).cos()); - let f_int = Vec3::uniform(-0.5 * (5.0f64).cos()); - let nabla_g = df_dt*0.1 + f*2.0 + f_int*0.4; - let actual_df_dt = solve_step_diff_eq( - nabla_g, - Vec3::uniform(0.1), - Vec3::uniform(2.0), - Vec3::uniform(0.4), - Vec3::zero(), - 0.0, // delta_t - f, // f_prev - f_int, // f_int_prev - Vec3::zero(), - ); - assert_vec3_eq(actual_df_dt, df_dt); - } - - // PORT-STATUS: skip - #[test] - fn step_diff_eq_3rd_order_no_dt() { - // nabla_g = 0.1 df/dt + 2 f + 0.4 f_int + 0.3 f_int_int - // Let f(t) = sin(2 t) - // then f'(t) = 2 cos(2 t) - // then \int f(t) = -1/2 cos(2 t) - // then \int \int f(t) = -1/4 sin(2 t) - let f = Vec3::unit()*(5.0f64).sin(); - let df_dt = Vec3::unit()*2.0 * (5.0f64).cos(); - let f_int = Vec3::unit()*-0.5 * (5.0f64).cos(); - let f_int_int = Vec3::unit()*-0.25 * (5.0f64).sin(); - let nabla_g = df_dt*0.1 + f*2.0 + f_int*0.4 + f_int_int*0.3; - let actual_df_dt = solve_step_diff_eq( - nabla_g, - Vec3::uniform(0.1), - Vec3::uniform(2.0), - Vec3::uniform(0.4), - Vec3::uniform(0.3), - 0.0, // delta_t - f, // f_prev - f_int, // f_int_prev - f_int_int, // f_int_int_prev - ); - assert_vec3_eq(actual_df_dt, df_dt); - } - - // PORT-STATUS: skip - #[test] - fn step_diff_eq_1st_order_homogenous() { - // 0 = f'(t) + f(t) - // solution: f(t) = k*e^-t - let f = |t: f32| -> Vec3 { - Vec3::unit()*(-t).exp() - }; - let df_dt = |t: f32| -> Vec3 { - Vec3::unit()*-(-t).exp() - }; - let f_int = |t: f32| -> Vec3 { - Vec3::unit()*-(-t).exp() - }; - let f_int_int = |t: f32| -> Vec3 { - Vec3::unit()*(-t).exp() - }; - let tl = 3.0000; - let tm = 3.0001; - let actual_df_dt = solve_step_diff_eq( - Vec3::zero(), - Vec3::unit(), // df_dt_coeff - Vec3::unit(), // f_coeff - Vec3::zero(), - Vec3::zero(), - tm - tl, // delta_t - f(tl), // f_prev - f_int(tl), - f_int_int(tl), - ); - assert_vec3_eq(actual_df_dt, df_dt(tm)); - } - - // PORT-STATUS: skip - #[test] - fn step_diff_eq_1st_order_nonhomogenous() { - // 1 = f'(t) + f(t) - // let f(t) = 1 + e^-t - let f = |t: f32| -> Vec3 { - Vec3::unit() + Vec3::unit()*(-t).exp() - }; - let df_dt = |t: f32| -> Vec3 { - Vec3::unit()*-(-t).exp() - }; - let f_int = |t: f32| -> Vec3 { - Vec3::unit()*-(-t).exp() - }; - let f_int_int = |t: f32| -> Vec3 { - Vec3::unit()*(-t).exp() - }; - let tl = 3.0000; - let tm = 3.0001; - let actual_df_dt = solve_step_diff_eq( - Vec3::unit(), // nabla_g - Vec3::unit(), // df_dt_coeff - Vec3::unit(), // f_coeff - Vec3::zero(), - Vec3::zero(), - tm - tl, // delta_t - f(tl), // f_prev - f_int(tl), - f_int_int(tl), - ); - assert_vec3_eq(actual_df_dt, df_dt(tm)); - } - - // PORT-STATUS: skip - #[test] - fn step_diff_eq_3rd_order_nonhomogenous() { - // c0 = af'(t) + bf(t) + c\int f(t) + d\int \int f(t) - // let f(t) = v0*e^wt, and \int \int f(t) absorbs C - // These constants are chosen arbitrarily: this test - // should pass for any "valid" values - let v0 = Vec3::new(3.0, 5.0, -0.5); - let w = Vec3::new(0.1, -2.0, 0.4); - let c0 = Vec3::new(-1.0, -2.0, -3.0); - let a = Vec3::new(0.1, 1.0, 2.0); // must be non-zero - let b = Vec3::new(0.5, -0.3, 0.0); - let c = Vec3::new(-5.0, 0.0, -0.1); - let d = -w.elem_mul(c + w.elem_mul(b + w.elem_mul(a))); - let f = |t: f64| -> Vec3 { - v0.elem_mul((w*t).exp()) - }; - let df_dt = |t: f64| -> Vec3 { - v0.elem_mul((w*t).exp()).elem_mul(w) - }; - let f_int = |t: f64| -> Vec3 { - v0.elem_mul((w*t).exp()).elem_div(w) - }; - let f_int_int = |t: f64| -> Vec3 { - v0.elem_mul((w*t).exp()).elem_div(w).elem_div(w) + c0.elem_div(d) - }; - let tl = 3.00000000; - let tm = 3.00000001; - let actual_df_dt = solve_step_diff_eq( - c0, // nabla_g - a, // df_dt_coeff - b, // f_coeff - c, // int_f_coeff - d, // int_int_f_coeff - tm - tl, // delta_t - f(tl), // f_prev - f_int(tl), // f_int_prev - f_int_int(tl), // f_int_int_prev - ); - assert_vec3_eq(actual_df_dt, df_dt(tm)); - } - - fn energy(s: &S) -> f32 { - 0.5 * s.feature_volume() * s.map_sum_enumerated(|_pos: Index, cell| { - // println!("{}: {}", _pos, cell.e()); - cell.e().mag_sq().to_f64() - }).to_f32() - } - - fn energy_now_and_then(state: &mut S, frames: u32) -> (f32, f32) { - let energy_0 = energy(state); - for _f in 0..frames { - // println!("e({}) = {}", _f, energy(state)); - state.step(); - } - let energy_1 = energy(state); - (energy_0, energy_1) - } - - // PORT-STATUS: done - #[test] - fn accessors() { - let mut state = SimState::>::new(Index((15, 17, 19).into()), 1e-6); - assert_eq!(state.width(), 15); - assert_eq!(state.height(), 17); - assert_eq!(state.depth(), 19); - assert_eq!(state.time(), 0.0); - // feat_size / C * courant: 1e-6 / 3e8 * 0.577 - assert_eq!(state.timestep(), 0.000000000000001924664829293337); - state.step(); - assert_eq!(state.time(), state.timestep()); - } - - struct ConstStim { - loc: Index, - e: Vec3, - h: Vec3, - } - impl Stimulus for ConstStim { - fn at(&self, _t_sec: R, _feat_size: R, loc: Index) -> Fields { - if loc == self.loc { - Fields { e: self.e.cast(), h: self.h.cast() } - } else { - Default::default() - } - } - } - - // PORT-STATUS: done - #[test] - fn energy_conservation_over_time() { - let mut state = SimState::>::new(Index((2001, 1, 1).into()), 1e-6); - for t in 0..100 { - let angle = (t as f32)/100.0*f32::two_pi(); - let sig = angle.sin(); - let gate = 0.5*(1.0 - angle.cos()); - //println!("{}", g.exp()); - state.step_multiple(1, &ConstStim { - loc: Index::new(0, 0, 1000), - e: Vec3::new_z(gate*sig/state.timestep()), - h: Vec3::zero(), - }); - } - let (energy_0, energy_1) = energy_now_and_then(&mut state, 800); - // This has some sensitivity to courant number. But surprisingly, not much to f32 v.s. f64 - assert_float_eq!(energy_1, energy_0, r2nd <= 1e-9); - } - - // PORT-STATUS: done - #[test] - fn sane_boundary_conditions() { - let mut state = SimState::>::new(Index((21, 21, 21).into()), 1e-6); - for t in 0..40 { - let angle = (t as f32)/10.0*f32::two_pi(); - let sig = angle.sin(); - let gate = 1e9 * 0.5*(1.0 - angle.cos()); - state.step_multiple(1, &ConstStim { - loc: Index::new(10, 10, 10), - e: Vec3::new_z(gate*sig/state.timestep()), - h: Vec3::zero(), - }); - } - let (energy_0, energy_1) = energy_now_and_then(&mut state, 500); - // Default boundary conditions reflect waves. - assert_float_eq!(energy_1, energy_0, r2nd <= 1e-2); - } - - /// Fill the world with the provided material and a stimulus. - /// Measure energy at the start, and then again after advancing many steps. - /// Return these two measurements (energy(t=0), energy(t=~=1000)) - fn conductor_test>>(mat: M) -> (f32, f32) { - let mut state = SimState::>::new(Index((201, 1, 1).into()), 1e-6); - state.fill_region(&WorldRegion, mat.into()); - for t in 0..100 { - let progress = (t as f32)/100.0; - let angle = 2.0*progress*f32::two_pi(); - let sig = angle.sin(); - let gate = (progress*f32::pi()).sin(); - state.step_multiple(1, &ConstStim { - loc: Index::new(100, 0, 0), - e: Vec3::new_z(gate*sig/state.timestep()), - h: Vec3::zero(), - }); - } - energy_now_and_then(&mut state, 1000) - } - - // PORT-STATUS: done - #[test] - fn conductor_dissipates_energy() { - let (energy_0, energy_1) = conductor_test(IsomorphicConductor::new(R64::from_f32(1e3))); - assert_float_eq!(energy_1/energy_0, 0.0, abs <= 1e-6); - } - - // PORT-STATUS: done - #[test] - fn anisotropic_conductor_inactive_x() { - let (energy_0, energy_1) = conductor_test(AnisomorphicConductor::new( - Vec3::new_x(1e3).cast() - )); - assert_gt!(energy_1, 0.9*energy_0); - // XXX: I think this gains energy because we only set the E field and not the H field - assert_lt!(energy_1, 1.5*energy_0); - } - - // PORT-STATUS: done - #[test] - fn anisotropic_conductor_inactive_y() { - let (energy_0, energy_1) = conductor_test(AnisomorphicConductor::new( - Vec3::new_y(1e3).cast() - )); - assert_gt!(energy_1, 0.9*energy_0); - assert_lt!(energy_1, 1.5*energy_0); - } - - // PORT-STATUS: done - #[test] - fn anisotropic_conductor_active_z() { - let (energy_0, energy_1) = conductor_test(AnisomorphicConductor::new( - Vec3::new_z(1e3).cast() - )); - assert_float_eq!(energy_1/energy_0, 0.0, abs <= 1e-6); - } - - fn state_for_pml(size: Index) -> SimState> { - let mut state = SimState::new(size, 1e-6); - let timestep = state.timestep(); - state.fill_boundary_using(size/4, |boundary_ness| { - let b = boundary_ness.elem_pow(3.0); - let conductivity = b * (0.5 / timestep); - // K scales up to 11; dc_shift scales from 0.05: - // https://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.461.4606&rep=rep1&type=pdf - // let propagation = Vec3::unit() + b * 10.0; - // let dc_shift = (Vec3::unit() - b) * 0.05; - // - // let propagation = Vec3::unit(); - // let dc_shift = Vec3::unit() * 1e-6; - Pml::new(conductivity) - }); - state - } - - /// Like state_for_pml, but the boundary is an ordinary electric conductor -- not PML - fn state_for_graded_conductor(size: Index) -> SimState> { - let mut state = SimState::new(size, 1e-6); - let timestep = state.timestep(); - state.fill_boundary_using(size/4, |boundary_ness| { - let b = boundary_ness.elem_pow(3.0); - let conductivity = f32::eps0() * b.mag() * 0.5 / timestep; - IsomorphicConductor::new(conductivity.cast()) - }); - state - } - - /// Stimulus which applies a (possibly different) e*sin(wt)*gate(t) at every point - struct PmlStim { - /// Maps location -> (field direction, frequency) - f: F, - t_end: f32, - sqrt_vol: f32, - } - impl (Vec3, f32) + Sync> Stimulus for PmlStim { - fn at(&self, t_sec: R, _feat_size: R, loc: Index) -> Fields { - let t_sec = t_sec.cast::(); - let angle = t_sec/self.t_end * f32::two_pi(); - let gate = 0.5*(1.0 - angle.cos()); - let (e, hz) = (self.f)(loc); - let sig_angle = angle*hz; - let sig = sig_angle.sin(); - // TODO: test H component? - Fields { - e: (e*(gate*sig / (self.t_end * self.sqrt_vol))).cast(), - h: Vec3::zero() - } - } - } - - /// 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)) - { - let stim = PmlStim { - f, - t_end: 50.0 * state.timestep(), - sqrt_vol: state.volume().sqrt(), - }; - for _t in 0..50 { - // println!("estim({}) = {}", _t, energy(state)); - state.step_multiple(1, &stim); - } - let (energy_0, energy_1) = energy_now_and_then(state, 1000); - // sanity check the starting energy - assert_gt!(energy_0, 1e-11); - assert_lt!(energy_0, 1.0); - energy_1 / energy_0 - } - - fn pml_test_over_region(state: &mut S, reg: R, f: F) -> f32 - where - R: Region, - F: Fn(Index) -> (Vec3, f32) + Sync, - S: AbstractSim, - { - let feat = state.feature_size(); - pml_test_full_interior(state, |idx| { - if reg.contains(idx.to_meters(feat)) { - f(idx) - } else { - (Vec3::zero(), 0.0) - } - }) - } - - fn pml_test_at(state: &mut S, e: Vec3, center: Index) -> f32 { - pml_test_full_interior(state, |idx| { - if idx == center { - (e, 1.0) - } else { - (Vec3::zero(), 1.0) - } - }) - } - - fn pml_test + Debug>(state: &mut S, e: Vec3, en_range: R) { - let center = Index(state.size().0 / 2); - let pml = pml_test_at(state, e, center); - // PML should absorb all energy - assert!(en_range.contains(&pml), "{} not in {:?}", pml, en_range); - } - - 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); - let en_baseline = pml_test_at(baseline_state, e, center); - // PML should absorb all energy - println!("en_pml: {}", en_pml); - println!("en_baseline: {}", en_baseline); - assert_float_eq!(en_pml, 0.0, abs <= 1.2e-6); - // XXX: why is PML energy sometimes MORE here? - assert_lt!(en_pml / en_baseline, 3.0, "{}/{}", en_pml, en_baseline); - // panic!("debugging"); - } - - // PORT-STATUS: skip - #[test] - fn pml_boundary_conditions_x() { - let mut pml_state = state_for_pml(Index::new(201, 1, 1)); - let mut baseline_state = state_for_graded_conductor(Index::new(201, 1, 1)); - pml_test_against_baseline(&mut pml_state, &mut baseline_state, Vec3::unit_z()); - } - - // PORT-STATUS: skip - #[test] - fn pml_boundary_conditions_y() { - let mut pml_state = state_for_pml(Index::new(1, 201, 1)); - let mut baseline_state = state_for_graded_conductor(Index::new(1, 201, 1)); - pml_test_against_baseline(&mut pml_state, &mut baseline_state, Vec3::unit_x()); - } - - // PORT-STATUS: skip - #[test] - fn pml_boundary_conditions_z() { - let mut pml_state = state_for_pml(Index::new(1, 1, 201)); - let mut baseline_state = state_for_graded_conductor(Index::new(1, 1, 201)); - pml_test_against_baseline(&mut pml_state, &mut baseline_state, Vec3::unit_y()); - } - - // PORT-STATUS: skip - #[test] - fn pml_boundary_conditions_plane_xy() { - let mut pml_state = state_for_pml(Index::new(201, 201, 1)); - let mut baseline_state = state_for_graded_conductor(Index::new(201, 201, 1)); - pml_test_against_baseline(&mut pml_state, &mut baseline_state, Vec3::unit_z()); - } - - // PORT-STATUS: skip - #[test] - fn pml_boundary_conditions_plane_xz() { - let mut pml_state = state_for_pml(Index::new(201, 1, 201)); - let mut baseline_state = state_for_graded_conductor(Index::new(201, 1, 201)); - pml_test_against_baseline(&mut pml_state, &mut baseline_state, Vec3::unit_y()); - } - - // PORT-STATUS: skip - #[test] - fn pml_boundary_conditions_plane_yz() { - let mut pml_state = state_for_pml(Index::new(1, 201, 201)); - let mut baseline_state = state_for_graded_conductor(Index::new(1, 201, 201)); - pml_test_against_baseline(&mut pml_state, &mut baseline_state, Vec3::unit_x()); - } - - fn state_for_monodirectional_pml(size: Index) -> SimState> { - let mut state = SimState::new(size, 1e-6); - let timestep = state.timestep(); - state.fill_boundary_using(size/4, |boundary_ness| { - let b = boundary_ness.elem_pow(3.0); - let coord_stretch = b * (0.5 / timestep); - // let coord_stretch = b * 0.5; - // Force the stretch to be only along one axis - let coord_stretch: Vec3<_> = match coord_stretch.to_tuple() { - (x, y, z) if x >= y && x >= z => (x, 0.0, 0.0), - (x, y, z) if y >= x && y >= z => (0.0, y, 0.0), - (x, y, z) if z >= x && z >= y => (0.0, 0.0, z), - _ => unreachable!(), - }.into(); - Pml::new(coord_stretch) - }); - state - } - - /// 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 - ) { - let mut state = SimState::>::new(size, 1e-6); - let timestep = state.timestep(); - state.fill_boundary_using(size/4, |boundary_ness| { - let b = boundary_ness.elem_pow(3.0); - 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); - Pml::::new(coord_stretch) - }); - let center = Index(state.size().0 / 2); - let en = pml_test_at(&mut state, e, center); - // XXX: would be better to compare it to an empty sim? - assert_float_eq!(en, 1.0, abs <= 0.1); - } - - // PORT-STATUS: skip - #[test] - fn pml_boundary_3d_inneffective_monodirectional_linear() { - for (size, e) in vec![ - (Index::new(201, 1, 1), Vec3::unit_y()), - (Index::new(201, 1, 1), Vec3::unit_z()), - (Index::new(1, 201, 1), Vec3::unit_x()), - (Index::new(1, 201, 1), Vec3::unit_z()), - (Index::new(1, 1, 201), Vec3::unit_x()), - (Index::new(1, 1, 201), Vec3::unit_y()), - ] { - pml_ineffective_mono_linear_test( - size, e, |v| Vec3::new(v.y(), v.z(), v.x()) - ); - pml_ineffective_mono_linear_test( - size, e, |v| Vec3::new(v.z(), v.x(), v.y()) - ); - } - } - - // PORT-STATUS: skip - #[test] - fn pml_boundary_3d_monodirectional_plane_xy() { - //let mut state = state_for_monodirectional_pml(Index::unit()*41); - let mut state = state_for_monodirectional_pml(Index::new(201, 201, 1)); - pml_test(&mut state, Vec3::unit_z(), ..7e-5); - } - - // PORT-STATUS: skip - #[test] - fn pml_boundary_3d_monodirectional_plane_xz() { - //let mut state = state_for_monodirectional_pml(Index::unit()*41); - let mut state = state_for_monodirectional_pml(Index::new(201, 1, 201)); - pml_test(&mut state, Vec3::unit_y(), ..7e-5); - } - - // PORT-STATUS: skip - #[test] - fn pml_boundary_3d_monodirectional_plane_yz() { - //let mut state = state_for_monodirectional_pml(Index::unit()*41); - let mut state = state_for_monodirectional_pml(Index::new(1, 201, 201)); - pml_test(&mut state, Vec3::unit_x(), ..7e-5); - } - - // PORT-STATUS: skip - #[test] - fn pml_boundary_3d_monodirectional_x() { - let mut state = state_for_monodirectional_pml(Index::unit()*41); - pml_test(&mut state, Vec3::unit_x(), 1e30..); - } - - // PORT-STATUS: skip - #[test] - fn pml_boundary_3d_monodirectional_y() { - let mut state = state_for_monodirectional_pml(Index::unit()*41); - pml_test(&mut state, Vec3::unit_y(), 1e30..); - } - - // PORT-STATUS: skip - #[test] - fn pml_boundary_3d_monodirectional_z() { - let mut state = state_for_monodirectional_pml(Index::unit()*41); - pml_test(&mut state, Vec3::unit_z(), 1e30..); - } - - // PORT-STATUS: skip - #[test] - fn pml_boundary_3d_multidirectional() { - for dir in vec![Vec3::unit_x(), Vec3::unit_y(), Vec3::unit_z()] { - let mut pml_state = state_for_pml(Index::unit()*41); - let mut baseline_state = state_for_graded_conductor(Index::unit()*41); - pml_test_against_baseline(&mut pml_state, &mut baseline_state, dir); - } - } - - // PORT-STATUS: skip - #[test] - fn pml_boundary_3d_multidirectional_non_axis_aligned() { - for dir in vec![ - Vec3::new(0.0, 1.0, 1.0), - Vec3::new(-1.0, 0.0, 1.0), - Vec3::new(1.0, -1.0, 0.0), - Vec3::new(1.0, -1.0, 1.0), - Vec3::new(-1.0, -1.0, -1.0), - Vec3::new(0.2, 0.5, 0.6), - Vec3::new(-2.0, 0.5, -0.5), - ] { - let mut pml_state = state_for_pml(Index::unit()*41); - let mut baseline_state = state_for_graded_conductor(Index::unit()*41); - pml_test_against_baseline(&mut pml_state, &mut baseline_state, dir); - } - } - - // PORT-STATUS: skip - #[test] - fn pml_boundary_3d_multidirectional_off_center() { - for dir in vec![Vec3::unit_x(), Vec3::unit_y(), Vec3::unit_z()] { - for center in vec![ - Index::new(15, 15, 15), - Index::new(15, 15, 25), - Index::new(15, 25, 15), - Index::new(15, 25, 25), - Index::new(25, 15, 15), - Index::new(25, 15, 25), - Index::new(25, 25, 15), - Index::new(25, 25, 25), - ] { - let mut pml_state = state_for_pml(Index::unit()*41); - let mut baseline_state = state_for_graded_conductor(Index::unit()*41); - let pml_en = pml_test_at(&mut pml_state, dir, center); - let baseline_en = pml_test_at(&mut baseline_state, dir, center); - assert_float_eq!(pml_en, 0.0, abs <= 2e-5); - // XXX: why is PML energy MORE here? - assert_lt!(pml_en / baseline_en, 2.0); - } - } - } - - // XXX the energies here are even worse than above. - // just don't use our PML implementation: it's not great. - // 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) - // }) - // } - // #[test] - // fn pml_boundary_3d_multidirectional_off_center_multisource() { - // for dir in vec![Vec3::unit_x(), Vec3::unit_y(), Vec3::unit_z()] { - // for cycles in vec![1.0, 2.0, 4.0, 8.0] { - // let mut pml_state = state_for_pml(Index::unit()*41); - // let mut baseline_state = state_for_graded_conductor(Index::unit()*41); - // let feat = pml_state.feature_size(); - // let interior = Cube::new(Index::new(15, 15, 15).to_meters(feat), Index::new(25, 25, 25).to_meters(feat)); - // let pml_en = pml_test_uniform_over_region(&mut pml_state, interior, dir, cycles); - // let baseline_en = pml_test_uniform_over_region(&mut baseline_state, interior, dir, cycles); - // assert_float_eq!(pml_en, 0.0, abs <= 4e-5); - // // XXX: why is PML energy MORE here? - // assert_lt!(pml_en / baseline_en, 5.0, "{} {} {}/{}", dir, cycles, pml_en, baseline_en); - // } - // } - // } - - // PORT-STATUS: skip - #[test] - fn pml_boundary_3d_chaotic_field() { - let excitation = |idx: Index| { - let seed = (idx.x() as u64) ^ ((idx.y() as u64) << 16) ^ ((idx.z() as u64) << 32); - let mut rng = StdRng::seed_from_u64(seed); - let dir = Vec3::new( - rng.gen_range(-1.0..1.0), - rng.gen_range(-1.0..1.0), - rng.gen_range(-1.0..1.0), - ); - // XXX only works if it's a whole number. I suppose this makes sense though, as - // other numbers would create higher harmonics when gated. - let hz = rng.gen_range(0..=5); - (dir, hz as _) - }; - let mut pml_state = state_for_pml(Index::unit()*41); - let mut baseline_state = state_for_graded_conductor(Index::unit()*41); - let feat = pml_state.feature_size(); - let interior = Cube::new(Index::new(15, 15, 15).to_meters(feat), Index::new(25, 25, 25).to_meters(feat)); - let en_pml = pml_test_over_region(&mut pml_state, interior, excitation); - let en_baseline = pml_test_over_region(&mut baseline_state, interior, excitation); - assert_float_eq!(en_pml, 0.0, abs <= 2e-5); - // XXX: why is PML energy MORE here? - assert_lt!(en_pml/en_baseline, 3.0); - } -} diff --git a/crates/coremem/src/sim/mod.rs b/crates/coremem/src/sim/mod.rs index 9b3d241..dbdd563 100644 --- a/crates/coremem/src/sim/mod.rs +++ b/crates/coremem/src/sim/mod.rs @@ -9,7 +9,6 @@ use rayon::prelude::*; use serde::{Serialize, Deserialize}; use std::iter::Sum; -pub mod legacy; pub mod spirv; pub mod units;