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

spirv_backend: use Real:: constants instead of inlined ones

Browse Source
This commit is contained in:
colin
2022-07-18 15:10:57 -07:00
parent 57338bcb4a
commit b8bcd68b98
3 changed files with 122 additions and 93 deletions
Download Patch File Download Diff File Expand all files Collapse all files

82
crates/spirv_backend/src/mat.rs
View File

@@ -1,6 +1,7 @@
use crate::support::Optional; use crate::support::Optional;
use coremem_types::mat::Material; use coremem_types::mat::Material;
use coremem_types::vec::Vec3; use coremem_types::vec::Vec3;
use coremem_types::real::{Real as _};
/// M(B) parallelogram /// M(B) parallelogram
@@ -74,23 +75,20 @@ pub struct MHPgram {
impl MHPgram { impl MHPgram {
/// h_intercept: X coordinate at which M is always zero. /// h_intercept: X coordinate at which M is always zero.
/// mu_r: relative mu value along the non-flat edges of the parallelogram. /// 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 { pub fn new(h_intercept: f32, mu_r: f32, max_m: f32) -> Self {
const MU0_INV: f32 = 795774.715025073;
let one_minus_mu_r_inv = 1.0 - 1.0/mu_r; let one_minus_mu_r_inv = 1.0 - 1.0/mu_r;
Self { 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, m_offset: h_intercept * one_minus_mu_r_inv,
max_m, max_m,
} }
} }
pub fn h_intercept(&self) -> f32 { pub fn h_intercept(&self) -> f32 {
const MU0_INV: f32 = 795774.715025073; f32::mu0_inv() * self.m_offset / self.b_mult
MU0_INV * self.m_offset / self.b_mult
} }
pub fn mu_r(&self) -> f32 { pub fn mu_r(&self) -> f32 {
const MU0: f32 = 1.2566370621219e-06; let mu_r_inv = 1.0 - f32::mu0()*self.b_mult;
let mu_r_inv = 1.0 - MU0*self.b_mult;
1.0 / mu_r_inv 1.0 / mu_r_inv
} }
@@ -177,8 +175,8 @@ pub struct Ferroxcube3R1MH;
impl Into<MHPgram> for Ferroxcube3R1MH { impl Into<MHPgram> for Ferroxcube3R1MH {
fn into(self) -> MHPgram { fn into(self) -> MHPgram {
const CURVE: MHPgram = MHPgram::new(25.0, 881.33, 44000.0); // TODO: how much (if any) penalty do we pay for this not being `const`?
CURVE MHPgram::new(25.0, 881.33, 44000.0)
} }
} }
@@ -312,55 +310,55 @@ mod test {
#[test] #[test]
fn mh_curve_edge_travel() { 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); 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(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(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(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, 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, 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(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(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(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(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(-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, -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(-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(-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(-300.0, 50.0*mu0), 0.0, 0.1); // interior (rightward) travel to H=50
} }
#[test] #[test]
fn mh_curve_interior() { fn mh_curve_interior() {
const MU0: f32 = 1.2566370621219e-06; let mu0 = f32::mu0();
let curve = MHPgram::new(50.0, 101.0, 500.0); let curve = MHPgram::new(50.0, 101.0, 500.0);
// max M (right) 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; // min M (right) happens at H=45; B = -455*mu0;
// max M (left) 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; // 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, 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, 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, -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, -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, -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, -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, -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, 5000.0*mu0), 500.0, 0.1); // rightward travel from H=-54 to H>>55
} }
#[test] #[test]
fn mh_curve_exterior() { fn mh_curve_exterior() {
const MU0: f32 = 1.2566370621219e-06; let mu0 = f32::mu0();
let curve = MHPgram::new(50.0, 101.0, 500.0); 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, 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 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(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(-501.0, 454.0*mu0), 400.0, 0.1); // exterior travel from H<<-55 (rounding error) to H=54
} }
} }

16
crates/spirv_backend/src/sim.rs
View File

@@ -3,6 +3,7 @@ use crate::support::{
Array3, Array3Mut, ArrayHandle, ArrayHandleMut, Optional, UnsizedArray Array3, Array3Mut, ArrayHandle, ArrayHandleMut, Optional, UnsizedArray
}; };
use coremem_types::mat::Material; use coremem_types::mat::Material;
use coremem_types::real::{Real as _};
use coremem_types::vec::{Vec3, Vec3u}; use coremem_types::vec::{Vec3, Vec3u};
#[derive(Copy, Clone)] #[derive(Copy, Clone)]
@@ -280,8 +281,7 @@ pub struct StepEContext<'a, M> {
impl<'a, M: Material<f32>> StepEContext<'a, M> { impl<'a, M: Material<f32>> StepEContext<'a, M> {
pub fn step_e(mut self) { pub fn step_e(mut self) {
// const EPS0_INV: f32 = 112940906737.1361; let twice_eps0 = f32::twice_eps0();
const TWICE_EPS0: f32 = 1.7708375625626e-11;
let deltas = self.in_h.delta_h(); let deltas = self.in_h.delta_h();
// \nabla x H // \nabla x H
let nabla_h = deltas.nabla() * self.inv_feature_size; let nabla_h = deltas.nabla() * self.inv_feature_size;
@@ -291,7 +291,7 @@ impl<'a, M: Material<f32>> StepEContext<'a, M> {
let sigma = self.mat.get_ref().conductivity(); let sigma = self.mat.get_ref().conductivity();
let e_prev = self.out_e.get(); let e_prev = self.out_e.get();
let delta_e = (nabla_h - e_prev.elem_mul(sigma)).elem_div( 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); )*(2.0*self.time_step);
// println!("spirv-step_e delta_e: {:?}", delta_e); // println!("spirv-step_e delta_e: {:?}", delta_e);
self.out_e.write(e_prev + delta_e + self.stim_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<f32>> StepHContext<'a, M> { impl<'a, M: Material<f32>> StepHContext<'a, M> {
pub fn step_h(mut self) { pub fn step_h(mut self) {
const MU0: f32 = 1.2566370621219e-06; let mu0 = f32::mu0();
const MU0_INV: f32 = 795774.715025073; let mu0_inv = f32::mu0_inv();
let deltas = self.in_e.delta_e(); let deltas = self.in_e.delta_e();
// println!("spirv-step_h delta_e_struct: {:?}", deltas); // println!("spirv-step_h delta_e_struct: {:?}", deltas);
// \nabla x E // \nabla x E
@@ -325,12 +325,12 @@ impl<'a, M: Material<f32>> StepHContext<'a, M> {
// Relation between these is: B = mu0*(H + M) // Relation between these is: B = mu0*(H + M)
let old_h = self.out_h.get(); let old_h = self.out_h.get();
let old_m = self.out_m.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 mat = self.mat.get_ref();
let new_m = mat.move_b_vec(old_m, new_b); 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); // println!("spirv-step_h delta_h: {:?}", delta_h);
self.out_h.write(new_h); self.out_h.write(new_h);
self.out_m.write(new_m); self.out_m.write(new_m);

117
crates/types/src/real.rs
View File

@@ -98,54 +98,82 @@ pub trait Real:
self == Self::zero() self == Self::zero()
} }
fn zero() -> Self { fn zero() -> Self;
Self::from_primitive(0) fn one() -> Self;
} fn two() -> Self;
fn tenth() -> Self { fn three() -> Self;
Self::one() / Self::ten() fn ten() -> Self;
} fn tenth() -> Self;
fn third() -> Self { fn third() -> Self;
Self::one() / Self::three() fn half() -> Self;
} fn pi() -> Self;
fn half() -> Self { fn two_pi() -> 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)
}
/// Speed of light in a vacuum; m/s. /// Speed of light in a vacuum; m/s.
/// Also equal to 1/sqrt(epsilon_0 mu_0) /// Also equal to 1/sqrt(epsilon_0 mu_0)
fn c() -> Self { fn c() -> Self;
Self::from_primitive(299792458)
}
/// Vaccum Permittivity /// Vaccum Permittivity
fn eps0() -> Self { fn eps0() -> Self;
Self::from_primitive(8.854187812813e-12) // F⋅m−1 fn twice_eps0() -> Self;
}
/// Vacuum Permeability /// Vacuum Permeability
fn mu0() -> Self { fn mu0() -> Self;
Self::from_primitive(1.2566370621219e-6) // H/m fn mu0_inv() -> Self;
} }
fn mu0_inv() -> Self {
Self::one() / Self::mu0() 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 { impl ToFloat for f32 {
@@ -157,6 +185,7 @@ impl ToFloat for f32 {
impl RealFeatures for f32 {} impl RealFeatures for f32 {}
impl Real for f32 { impl Real for f32 {
decl_consts!(Real::from_f32);
fn from_primitive<P: ToFloat>(p: P) -> Self { fn from_primitive<P: ToFloat>(p: P) -> Self {
Self::from_f32(p.to_f32()) Self::from_f32(p.to_f32())
} }
@@ -201,6 +230,7 @@ impl ToFloat for f64 {
impl RealFeatures for f64 {} impl RealFeatures for f64 {}
impl Real for f64 { impl Real for f64 {
decl_consts!(Real::from_f64);
fn from_f64(f: f64) -> Self { fn from_f64(f: f64) -> Self {
f f
} }
@@ -347,6 +377,7 @@ impl<T: ToFloat> ToFloat for Finite<T> {
impl<T: RealFeatures> RealFeatures for Finite<T> {} impl<T: RealFeatures> RealFeatures for Finite<T> {}
impl<T: Real> Real for Finite<T> { impl<T: Real> Real for Finite<T> {
decl_consts!(Self::from_primitive);
fn from_primitive<P: ToFloat>(p: P) -> Self { fn from_primitive<P: ToFloat>(p: P) -> Self {
Self::new(T::from_primitive(p)) Self::new(T::from_primitive(p))
} }