move MBPgram, MHPgram out of spirv_backend into coremem_types::mat
later this can be shared with CPU backend.
This commit is contained in:
@@ -1,157 +1,13 @@
|
||||
use crate::support::Optional;
|
||||
use coremem_types::mat::Material;
|
||||
use coremem_types::mat::{Material, MBPgram, MHPgram};
|
||||
use coremem_types::vec::Vec3;
|
||||
use coremem_types::real::{Real as _};
|
||||
|
||||
|
||||
/// M(B) parallelogram
|
||||
///
|
||||
///```ignore
|
||||
/// ____________
|
||||
/// / /
|
||||
/// / . /
|
||||
/// / /
|
||||
/// /___________/
|
||||
/// ```
|
||||
///
|
||||
/// The `.` depicts (0, 0). X axis is B; y axis is M.
|
||||
/// As B increases, M remains constant until it hits an edge.
|
||||
/// Then M rises up to its max.
|
||||
/// Same thing happens on the left edge, as B decreases and M falls to its min.
|
||||
#[derive(Copy, Clone, Default, PartialEq)]
|
||||
pub struct MBPgram {
|
||||
/// X coordinate at which the upward slope starts
|
||||
pub b_start: f32,
|
||||
/// X coordinate at which the upward slope ends
|
||||
pub b_end: f32,
|
||||
/// Vertical range of the graph
|
||||
pub max_m: f32,
|
||||
}
|
||||
|
||||
/// f(x0) = -ymax
|
||||
/// f(x1) = ymax
|
||||
fn eval_bounded_line(x: f32, x0: f32, x1: f32, ymax: f32) -> f32 {
|
||||
let x = x - x0;
|
||||
let x1 = x1 - x0;
|
||||
// Now, f(0) = -ymax, f(x1) = ymax
|
||||
let x = x.max(0.0).min(x1);
|
||||
ymax * (2.0*x/x1 - 1.0)
|
||||
}
|
||||
|
||||
impl MBPgram {
|
||||
pub fn new(b_start: f32, b_end: f32, max_m: f32) -> Self {
|
||||
Self { b_start, b_end, max_m }
|
||||
}
|
||||
|
||||
/// Return the new `M`
|
||||
pub fn move_b(self, m: f32, target_b: f32) -> f32 {
|
||||
// Evaluate the curve on the right:
|
||||
let y_right = eval_bounded_line(target_b, self.b_start, self.b_end, self.max_m);
|
||||
let y_left = eval_bounded_line(target_b, -self.b_end, -self.b_start, self.max_m);
|
||||
m.max(y_right).min(y_left)
|
||||
}
|
||||
}
|
||||
|
||||
impl Material<f32> for MBPgram {
|
||||
fn move_b_vec(&self, m: Vec3<f32>, target_b: Vec3<f32>) -> Vec3<f32> {
|
||||
Vec3::new(
|
||||
self.move_b(m.x(), target_b.x()),
|
||||
self.move_b(m.y(), target_b.y()),
|
||||
self.move_b(m.z(), target_b.z()),
|
||||
)
|
||||
}
|
||||
}
|
||||
|
||||
#[derive(Copy, Clone, Debug, Default, PartialEq)]
|
||||
pub struct MHPgram {
|
||||
/// optimized form of mu_0^-1 * (1-mu_r^-1)
|
||||
b_mult: f32,
|
||||
/// optimized form of h_intercept * (1-mu_r^-1)
|
||||
m_offset: f32,
|
||||
/// Vertical range of the graph
|
||||
pub max_m: f32,
|
||||
}
|
||||
|
||||
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 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: 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 {
|
||||
f32::mu0_inv() * self.m_offset / self.b_mult
|
||||
}
|
||||
pub fn mu_r(&self) -> f32 {
|
||||
let mu_r_inv = 1.0 - f32::mu0()*self.b_mult;
|
||||
1.0 / mu_r_inv
|
||||
}
|
||||
|
||||
/// Return the new `M`
|
||||
pub fn move_b(self, m: f32, target_b: f32) -> f32 {
|
||||
// The point may exist outside the parallelogram.
|
||||
// The right slope is defined by:
|
||||
// B(H)/mu0 = h_intercept + (h-h_intercept)*mu_r
|
||||
// "unoptimized" solution:
|
||||
// let target_mh = target_b*MU0_INV;
|
||||
// let h_on_right = (target_mh-self.h_intercept)/self.mu_r + self.h_intercept;
|
||||
// let m_on_right = target_mh - h_on_right;
|
||||
// Left:
|
||||
// B(H)/mu0 = -h_intercept + (h+h_intercept)*mu_r
|
||||
// "unoptimized" solution:
|
||||
// let h_on_left = (target_mh+self.h_intercept)/self.mu_r - self.h_intercept;
|
||||
// let m_on_left = target_mh - h_on_left;
|
||||
|
||||
let m_on_right = target_b*self.b_mult - self.m_offset;
|
||||
let m_on_left = target_b*self.b_mult + self.m_offset;
|
||||
|
||||
if m_on_right >= m {
|
||||
// if m_on_right <= self.max_m {
|
||||
// // rightward edge movement
|
||||
// m_on_right
|
||||
// } else {
|
||||
// // right of saturation
|
||||
// self.max_m
|
||||
// }
|
||||
m_on_right.min(self.max_m)
|
||||
} else if m_on_left <= m {
|
||||
// if m_on_left >= -self.max_m {
|
||||
// // leftward edge movement
|
||||
// m_on_left
|
||||
// } else {
|
||||
// // left of saturation
|
||||
// -self.max_m
|
||||
// }
|
||||
m_on_left.max(-self.max_m)
|
||||
} else {
|
||||
// interior movement
|
||||
m
|
||||
}
|
||||
// this boils down to:
|
||||
// m.clamp(m_on_left, m_on_right).clamp(-max_m, max_m)
|
||||
}
|
||||
}
|
||||
|
||||
impl Material<f32> for MHPgram {
|
||||
fn move_b_vec(&self, m: Vec3<f32>, target_b: Vec3<f32>) -> Vec3<f32> {
|
||||
Vec3::new(
|
||||
self.move_b(m.x(), target_b.x()),
|
||||
self.move_b(m.y(), target_b.y()),
|
||||
self.move_b(m.z(), target_b.z()),
|
||||
)
|
||||
}
|
||||
}
|
||||
|
||||
#[derive(Copy, Clone, Default, PartialEq)]
|
||||
pub struct FullyGenericMaterial {
|
||||
pub conductivity: Vec3<f32>,
|
||||
pub m_b_curve: Optional<MBPgram>,
|
||||
pub m_h_curve: Optional<MHPgram>,
|
||||
pub m_b_curve: Optional<MBPgram<f32>>,
|
||||
pub m_h_curve: Optional<MHPgram<f32>>,
|
||||
}
|
||||
|
||||
impl Material<f32> for FullyGenericMaterial {
|
||||
@@ -173,8 +29,8 @@ impl Material<f32> for FullyGenericMaterial {
|
||||
#[derive(Copy, Clone, Default, PartialEq)]
|
||||
pub struct Ferroxcube3R1MH;
|
||||
|
||||
impl Into<MHPgram> for Ferroxcube3R1MH {
|
||||
fn into(self) -> MHPgram {
|
||||
impl Into<MHPgram<f32>> for Ferroxcube3R1MH {
|
||||
fn into(self) -> MHPgram<f32> {
|
||||
// TODO: how much (if any) penalty do we pay for this not being `const`?
|
||||
MHPgram::new(25.0, 881.33, 44000.0)
|
||||
}
|
||||
@@ -182,7 +38,7 @@ impl Into<MHPgram> for Ferroxcube3R1MH {
|
||||
|
||||
impl Material<f32> for Ferroxcube3R1MH {
|
||||
fn move_b_vec(&self, m: Vec3<f32>, target_b: Vec3<f32>) -> Vec3<f32> {
|
||||
let curve: MHPgram = (*self).into();
|
||||
let curve: MHPgram<f32> = (*self).into();
|
||||
curve.move_b_vec(m, target_b)
|
||||
}
|
||||
}
|
||||
@@ -225,140 +81,3 @@ impl<M: Material<f32> + Copy> Material<f32> for IsoConductorOr<M> {
|
||||
}
|
||||
}
|
||||
|
||||
#[cfg(test)]
|
||||
mod test {
|
||||
use super::*;
|
||||
|
||||
fn assert_eq_approx(actual: f32, expected: f32, max_error: f32) {
|
||||
assert!((actual - expected).abs() < max_error, "{} != {}", actual, expected);
|
||||
}
|
||||
|
||||
fn assert_eq_approx_tight(actual: f32, expected: f32) {
|
||||
assert_eq_approx(actual, expected, 1e-5);
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn mb_curve_interior() {
|
||||
let curve = MBPgram::new(4.0, 6.0, 20.0f32);
|
||||
|
||||
assert_eq_approx_tight(curve.move_b(0.0, 2.0), 0.0);
|
||||
assert_eq_approx_tight(curve.move_b(0.0, 5.0), 0.0);
|
||||
assert_eq_approx_tight(curve.move_b(1.0, 5.0), 1.0);
|
||||
assert_eq_approx_tight(curve.move_b(20.0, 5.0), 20.0);
|
||||
assert_eq_approx_tight(curve.move_b(-20.0, 4.0), -20.0);
|
||||
assert_eq_approx_tight(curve.move_b(-20.0, -6.0), -20.0);
|
||||
assert_eq_approx_tight(curve.move_b(20.0, -4.0), 20.0);
|
||||
assert_eq_approx_tight(curve.move_b(10.0, -2.0), 10.0);
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn mb_curve_exterior() {
|
||||
let curve = MBPgram::new(4.0, 6.0, 20.0f32);
|
||||
assert_eq_approx_tight(curve.move_b(0.0, 6.0), 20.0);
|
||||
assert_eq_approx_tight(curve.move_b(0.0, 7.0), 20.0);
|
||||
assert_eq_approx_tight(curve.move_b(0.0, -6.0), -20.0);
|
||||
assert_eq_approx_tight(curve.move_b(0.0, -7.0), -20.0);
|
||||
|
||||
assert_eq_approx_tight(curve.move_b(2.0, -6.0), -20.0);
|
||||
assert_eq_approx_tight(curve.move_b(20.0, -6.0), -20.0);
|
||||
assert_eq_approx_tight(curve.move_b(20.0, -5.0), 0.0);
|
||||
assert_eq_approx_tight(curve.move_b(20.0, -4.5), 10.0);
|
||||
assert_eq_approx_tight(curve.move_b(-15.0, 4.5), -10.0);
|
||||
assert_eq_approx_tight(curve.move_b(-15.0, 5.0), 0.0);
|
||||
assert_eq_approx_tight(curve.move_b(-15.0, 5.5), 10.0);
|
||||
assert_eq_approx_tight(curve.move_b(-15.0, 7.5), 20.0);
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn mb_curve_3r1() {
|
||||
// slope of 3r1 is about M=793210*B
|
||||
// This is almost identical to iron (795615)!
|
||||
let curve = MBPgram::new(-0.3899, 0.3900, 310000f32);
|
||||
|
||||
// magnetizing:
|
||||
// v.s. 198893 in B(H) curve
|
||||
assert_eq_approx(curve.move_b(0.0, 0.250), 198703.0, 1.0);
|
||||
// v.s. 278321 in B(H) curve
|
||||
assert_eq_approx(curve.move_b(198703.0, 0.350), 278201.0, 1.0);
|
||||
assert_eq_approx(curve.move_b(278201.0, 0.390), 310000.0, 1.0);
|
||||
|
||||
// de-magnetizing:
|
||||
// From saturation, decreasing B causes NO decrease in M: instead, it causes a decrease in
|
||||
// H. This is probably BAD: in the B(H) curve, a large change in H always causes a large
|
||||
// change in B. The movement of H here is likely to induce current, whereas it SHOULDN'T.
|
||||
assert_eq_approx(curve.move_b(310000.0, 0.38995), 310000.0, 1.0);
|
||||
// v.s. 258626 in B(H); H = 220
|
||||
assert_eq_approx(curve.move_b(310000.0, 0.325), 258406.0, 1.0);
|
||||
// here's where H crosses 0 (v.s. B=0.325 in the B(H) curve... quite a difference)
|
||||
assert_eq_approx(curve.move_b(310000.0, 0.050), 39788.438, 1.0);
|
||||
// v.s. 35.0 in B(H)
|
||||
assert_eq_approx(curve.move_b(310000.0, 0.0), 39.75, 1.0);
|
||||
|
||||
// negative magnetization:
|
||||
// erase the magnetization: H = -40
|
||||
assert_eq_approx(curve.move_b(310000.0, -0.00005), 0.0, 1.0);
|
||||
// the magnetization has been completely erased:
|
||||
assert_eq_approx(curve.move_b(310000.0, -0.25), -198703.0, 1.0);
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn mh_curve_parameters() {
|
||||
let curve = MHPgram::new(50.0, 101.0, 500.0);
|
||||
assert_eq_approx(curve.h_intercept(), 50.0, 1e-3);
|
||||
assert_eq_approx(curve.mu_r(), 101.0, 1e-3);
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn mh_curve_edge_travel() {
|
||||
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
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn mh_curve_interior() {
|
||||
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;
|
||||
|
||||
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
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn mh_curve_exterior() {
|
||||
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, -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
|
||||
|
||||
}
|
||||
}
|
||||
|
||||
Reference in New Issue
Block a user