diff --git a/crates/spirv_backend/src/mat.rs b/crates/spirv_backend/src/mat.rs index 0aa3c81..e4cab64 100644 --- a/crates/spirv_backend/src/mat.rs +++ b/crates/spirv_backend/src/mat.rs @@ -1,6 +1,7 @@ use crate::support::Optional; use coremem_types::mat::Material; use coremem_types::vec::Vec3; +use coremem_types::real::{Real as _}; /// M(B) parallelogram @@ -74,23 +75,20 @@ pub struct MHPgram { impl MHPgram { /// h_intercept: X coordinate at which M is always zero. /// mu_r: relative mu value along the non-flat edges of the parallelogram. - pub const fn new(h_intercept: f32, mu_r: f32, max_m: f32) -> Self { - const MU0_INV: f32 = 795774.715025073; + pub fn new(h_intercept: f32, mu_r: f32, max_m: f32) -> Self { let one_minus_mu_r_inv = 1.0 - 1.0/mu_r; Self { - b_mult: MU0_INV * one_minus_mu_r_inv, + b_mult: f32::mu0_inv() * one_minus_mu_r_inv, m_offset: h_intercept * one_minus_mu_r_inv, max_m, } } pub fn h_intercept(&self) -> f32 { - const MU0_INV: f32 = 795774.715025073; - MU0_INV * self.m_offset / self.b_mult + f32::mu0_inv() * self.m_offset / self.b_mult } pub fn mu_r(&self) -> f32 { - const MU0: f32 = 1.2566370621219e-06; - let mu_r_inv = 1.0 - MU0*self.b_mult; + let mu_r_inv = 1.0 - f32::mu0()*self.b_mult; 1.0 / mu_r_inv } @@ -177,8 +175,8 @@ pub struct Ferroxcube3R1MH; impl Into for Ferroxcube3R1MH { fn into(self) -> MHPgram { - const CURVE: MHPgram = MHPgram::new(25.0, 881.33, 44000.0); - CURVE + // TODO: how much (if any) penalty do we pay for this not being `const`? + MHPgram::new(25.0, 881.33, 44000.0) } } @@ -312,55 +310,55 @@ mod test { #[test] fn mh_curve_edge_travel() { - const MU0: f32 = 1.2566370621219e-06; + let mu0 = f32::mu0(); let curve = MHPgram::new(50.0, 101.0, 500.0); - assert_eq_approx(curve.move_b(0.0, 151.0*MU0), 100.0, 0.1); // rightward travel along edge - assert_eq_approx(curve.move_b(100.0, 252.0*MU0), 200.0, 0.1); - assert_eq_approx(curve.move_b(200.0, 555.0*MU0), 500.0, 0.1); - assert_eq_approx(curve.move_b(500.0, 500.0*MU0), 500.0, 0.1); // back (leftward) travel to H=0 - assert_eq_approx(curve.move_b(500.0, 455.0*MU0), 500.0, 0.1); // back (leftward) travel to H=-45 - assert_eq_approx(curve.move_b(500.0, 354.0*MU0), 400.0, 0.1); // back (leftward) travel to H=-46 - assert_eq_approx(curve.move_b(400.0, 253.0*MU0), 300.0, 0.1); // back (leftward) travel to H=-47 - assert_eq_approx(curve.move_b(300.0, -50.0*MU0), 0.0, 0.1); // back (leftward) travel to H=-50 - assert_eq_approx(curve.move_b(0.0, -151.0*MU0), -100.0, 0.1); // back (leftward) travel to H=-51 - assert_eq_approx(curve.move_b(-100.0, -555.0*MU0), -500.0, 0.1); // back (leftward) travel to H=-55 - assert_eq_approx(curve.move_b(-500.0, -456.0*MU0), -500.0, 0.1); // interior (rightward) travel to H=44 - assert_eq_approx(curve.move_b(-500.0, -354.0*MU0), -400.0, 0.1); // interior (rightward) travel to H=46 - assert_eq_approx(curve.move_b(-400.0, -253.0*MU0), -300.0, 0.1); // interior (rightward) travel to H=47 - assert_eq_approx(curve.move_b(-300.0, 50.0*MU0), 0.0, 0.1); // interior (rightward) travel to H=50 + assert_eq_approx(curve.move_b(0.0, 151.0*mu0), 100.0, 0.1); // rightward travel along edge + assert_eq_approx(curve.move_b(100.0, 252.0*mu0), 200.0, 0.1); + assert_eq_approx(curve.move_b(200.0, 555.0*mu0), 500.0, 0.1); + assert_eq_approx(curve.move_b(500.0, 500.0*mu0), 500.0, 0.1); // back (leftward) travel to H=0 + assert_eq_approx(curve.move_b(500.0, 455.0*mu0), 500.0, 0.1); // back (leftward) travel to H=-45 + assert_eq_approx(curve.move_b(500.0, 354.0*mu0), 400.0, 0.1); // back (leftward) travel to H=-46 + assert_eq_approx(curve.move_b(400.0, 253.0*mu0), 300.0, 0.1); // back (leftward) travel to H=-47 + assert_eq_approx(curve.move_b(300.0, -50.0*mu0), 0.0, 0.1); // back (leftward) travel to H=-50 + assert_eq_approx(curve.move_b(0.0, -151.0*mu0), -100.0, 0.1); // back (leftward) travel to H=-51 + assert_eq_approx(curve.move_b(-100.0, -555.0*mu0), -500.0, 0.1); // back (leftward) travel to H=-55 + assert_eq_approx(curve.move_b(-500.0, -456.0*mu0), -500.0, 0.1); // interior (rightward) travel to H=44 + assert_eq_approx(curve.move_b(-500.0, -354.0*mu0), -400.0, 0.1); // interior (rightward) travel to H=46 + assert_eq_approx(curve.move_b(-400.0, -253.0*mu0), -300.0, 0.1); // interior (rightward) travel to H=47 + assert_eq_approx(curve.move_b(-300.0, 50.0*mu0), 0.0, 0.1); // interior (rightward) travel to H=50 } #[test] fn mh_curve_interior() { - const MU0: f32 = 1.2566370621219e-06; + let mu0 = f32::mu0(); let curve = MHPgram::new(50.0, 101.0, 500.0); - // max M (right) happens at H=55; B = 555*MU0; - // min M (right) happens at H=45; B = -455*MU0; - // max M (left) happens at H=-45; B = 455*MU0; - // min M (left) happens at H=-55; B = -555*MU0; + // max M (right) happens at H=55; B = 555*mu0; + // min M (right) happens at H=45; B = -455*mu0; + // max M (left) happens at H=-45; B = 455*mu0; + // min M (left) happens at H=-55; B = -555*mu0; - assert_eq_approx(curve.move_b(0.0, 40.0*MU0), 0.0, 0.1); - assert_eq_approx(curve.move_b(0.0, 50.0*MU0), 0.0, 0.1); - assert_eq_approx(curve.move_b(0.0, -40.0*MU0), 0.0, 0.1); - assert_eq_approx(curve.move_b(0.0, -50.0*MU0), 0.0, 0.1); + assert_eq_approx(curve.move_b(0.0, 40.0*mu0), 0.0, 0.1); + assert_eq_approx(curve.move_b(0.0, 50.0*mu0), 0.0, 0.1); + assert_eq_approx(curve.move_b(0.0, -40.0*mu0), 0.0, 0.1); + assert_eq_approx(curve.move_b(0.0, -50.0*mu0), 0.0, 0.1); - assert_eq_approx(curve.move_b(-400.0, -400.0*MU0), -400.0, 0.1); // rightward travel from H=-54 to NOT H=-53.5 - assert_eq_approx(curve.move_b(-400.0, -355.0*MU0), -400.0, 0.1); // rightward travel from H=-54 to H=45 - assert_eq_approx(curve.move_b(-400.0, -354.0*MU0), -400.0, 0.1); // rightward travel from H=-54 to H=46 - assert_eq_approx(curve.move_b(-400.0, 5000.0*MU0), 500.0, 0.1); // rightward travel from H=-54 to H>>55 + assert_eq_approx(curve.move_b(-400.0, -400.0*mu0), -400.0, 0.1); // rightward travel from H=-54 to NOT H=-53.5 + assert_eq_approx(curve.move_b(-400.0, -355.0*mu0), -400.0, 0.1); // rightward travel from H=-54 to H=45 + assert_eq_approx(curve.move_b(-400.0, -354.0*mu0), -400.0, 0.1); // rightward travel from H=-54 to H=46 + assert_eq_approx(curve.move_b(-400.0, 5000.0*mu0), 500.0, 0.1); // rightward travel from H=-54 to H>>55 } #[test] fn mh_curve_exterior() { - const MU0: f32 = 1.2566370621219e-06; + let mu0 = f32::mu0(); let curve = MHPgram::new(50.0, 101.0, 500.0); - assert_eq_approx(curve.move_b(500.0, 556.0*MU0), 500.0, 0.1); // exterior travel to H=56 - assert_eq_approx(curve.move_b(500.0, 5000.0*MU0), 500.0, 0.1); // exterior travel to H>>55 + assert_eq_approx(curve.move_b(500.0, 556.0*mu0), 500.0, 0.1); // exterior travel to H=56 + assert_eq_approx(curve.move_b(500.0, 5000.0*mu0), 500.0, 0.1); // exterior travel to H>>55 - assert_eq_approx(curve.move_b(500.0, -5000.0*MU0), -500.0, 0.1); // exterior travel from H>>55 to H << -55 - assert_eq_approx(curve.move_b(-501.0, 454.0*MU0), 400.0, 0.1); // exterior travel from H<<-55 (rounding error) to H=54 + assert_eq_approx(curve.move_b(500.0, -5000.0*mu0), -500.0, 0.1); // exterior travel from H>>55 to H << -55 + assert_eq_approx(curve.move_b(-501.0, 454.0*mu0), 400.0, 0.1); // exterior travel from H<<-55 (rounding error) to H=54 } } diff --git a/crates/spirv_backend/src/sim.rs b/crates/spirv_backend/src/sim.rs index 17c9a28..6158d47 100644 --- a/crates/spirv_backend/src/sim.rs +++ b/crates/spirv_backend/src/sim.rs @@ -3,6 +3,7 @@ use crate::support::{ Array3, Array3Mut, ArrayHandle, ArrayHandleMut, Optional, UnsizedArray }; use coremem_types::mat::Material; +use coremem_types::real::{Real as _}; use coremem_types::vec::{Vec3, Vec3u}; #[derive(Copy, Clone)] @@ -280,8 +281,7 @@ pub struct StepEContext<'a, M> { impl<'a, M: Material> StepEContext<'a, M> { pub fn step_e(mut self) { - // const EPS0_INV: f32 = 112940906737.1361; - const TWICE_EPS0: f32 = 1.7708375625626e-11; + let twice_eps0 = f32::twice_eps0(); let deltas = self.in_h.delta_h(); // \nabla x H let nabla_h = deltas.nabla() * self.inv_feature_size; @@ -291,7 +291,7 @@ impl<'a, M: Material> StepEContext<'a, M> { let sigma = self.mat.get_ref().conductivity(); let e_prev = self.out_e.get(); let delta_e = (nabla_h - e_prev.elem_mul(sigma)).elem_div( - sigma*self.time_step + Vec3::uniform(TWICE_EPS0) + sigma*self.time_step + Vec3::uniform(twice_eps0) )*(2.0*self.time_step); // println!("spirv-step_e delta_e: {:?}", delta_e); self.out_e.write(e_prev + delta_e + self.stim_e); @@ -313,8 +313,8 @@ pub struct StepHContext<'a, M> { impl<'a, M: Material> StepHContext<'a, M> { pub fn step_h(mut self) { - const MU0: f32 = 1.2566370621219e-06; - const MU0_INV: f32 = 795774.715025073; + let mu0 = f32::mu0(); + let mu0_inv = f32::mu0_inv(); let deltas = self.in_e.delta_e(); // println!("spirv-step_h delta_e_struct: {:?}", deltas); // \nabla x E @@ -325,12 +325,12 @@ impl<'a, M: Material> StepHContext<'a, M> { // Relation between these is: B = mu0*(H + M) let old_h = self.out_h.get(); let old_m = self.out_m.get(); - let old_b = (old_h + old_m) * MU0; + let old_b = (old_h + old_m) * mu0; - let new_b = old_b + delta_b + self.stim_h * MU0; + let new_b = old_b + delta_b + self.stim_h * mu0; let mat = self.mat.get_ref(); let new_m = mat.move_b_vec(old_m, new_b); - let new_h = new_b * MU0_INV - new_m; + let new_h = new_b * mu0_inv - new_m; // println!("spirv-step_h delta_h: {:?}", delta_h); self.out_h.write(new_h); self.out_m.write(new_m); diff --git a/crates/types/src/real.rs b/crates/types/src/real.rs index c21d2df..8c6ba7e 100644 --- a/crates/types/src/real.rs +++ b/crates/types/src/real.rs @@ -98,54 +98,82 @@ pub trait Real: self == Self::zero() } - fn zero() -> Self { - Self::from_primitive(0) - } - fn tenth() -> Self { - Self::one() / Self::ten() - } - fn third() -> Self { - Self::one() / Self::three() - } - fn half() -> Self { - Self::from_f32(0.5) - } - fn one() -> Self { - Self::from_primitive(1) - } - fn two() -> Self { - Self::from_f32(2.0) - } - fn three() -> Self { - Self::from_f32(3.0) - } - fn ten() -> Self { - Self::from_f32(10.0) - } - - fn pi() -> Self { - Self::from_f64(core::f64::consts::PI) - } - fn two_pi() -> Self { - Self::from_f64(core::f64::consts::PI * 2.0) - } + fn zero() -> Self; + fn one() -> Self; + fn two() -> Self; + fn three() -> Self; + fn ten() -> Self; + fn tenth() -> Self; + fn third() -> Self; + fn half() -> Self; + fn pi() -> Self; + fn two_pi() -> Self; /// Speed of light in a vacuum; m/s. /// Also equal to 1/sqrt(epsilon_0 mu_0) - fn c() -> Self { - Self::from_primitive(299792458) - } + fn c() -> Self; /// Vaccum Permittivity - fn eps0() -> Self { - Self::from_primitive(8.854187812813e-12) // F⋅m−1 - } + fn eps0() -> Self; + fn twice_eps0() -> Self; /// Vacuum Permeability - fn mu0() -> Self { - Self::from_primitive(1.2566370621219e-6) // H/m - } - fn mu0_inv() -> Self { - Self::one() / Self::mu0() - } + fn mu0() -> Self; + fn mu0_inv() -> Self; +} + +macro_rules! decl_consts { + ($wrap:expr) => { + fn zero() -> Self { + $wrap(0.0) + } + fn one() -> Self { + $wrap(1.0) + } + fn two() -> Self { + $wrap(2.0) + } + fn three() -> Self { + $wrap(3.0) + } + fn ten() -> Self { + $wrap(10.0) + } + fn tenth() -> Self { + $wrap(0.1) + } + fn third() -> Self { + $wrap(0.3333333333333333) + } + fn half() -> Self { + $wrap(0.5) + } + + fn pi() -> Self { + $wrap(3.141592653589793) + } + fn two_pi() -> Self { + $wrap(6.283185307179586) + } + + /// Speed of light in a vacuum; m/s. + /// Also equal to 1/sqrt(epsilon_0 mu_0) + fn c() -> Self { + $wrap(299792458.0) + } + /// Vaccum Permittivity + fn eps0() -> Self { + $wrap(8.854187812813e-12) // F⋅m−1 + } + fn twice_eps0() -> Self { + $wrap(1.7708375625626e-11) // F⋅m−1 + } + /// Vacuum Permeability + fn mu0() -> Self { + $wrap(1.2566370621219e-6) // H/m + } + fn mu0_inv() -> Self { + $wrap(795774.715025073) // m/H + } + }; } impl ToFloat for f32 { @@ -157,6 +185,7 @@ impl ToFloat for f32 { impl RealFeatures for f32 {} impl Real for f32 { + decl_consts!(Real::from_f32); fn from_primitive(p: P) -> Self { Self::from_f32(p.to_f32()) } @@ -201,6 +230,7 @@ impl ToFloat for f64 { impl RealFeatures for f64 {} impl Real for f64 { + decl_consts!(Real::from_f64); fn from_f64(f: f64) -> Self { f } @@ -347,6 +377,7 @@ impl ToFloat for Finite { impl RealFeatures for Finite {} impl Real for Finite { + decl_consts!(Self::from_primitive); fn from_primitive(p: P) -> Self { Self::new(T::from_primitive(p)) }