fdtd-coremem/crates/cross/src/mat/mb_pgram.rs

170 lines
6.1 KiB
Rust

use crate::compound::enumerated::{Discr, DiscriminantCodable};
use crate::compound::peano::Peano;
use crate::mat::Material;
use crate::real::Real;
use crate::vec::Vec3;
#[cfg(feature = "serde")]
use serde::{Serialize, Deserialize};
/// M(B) parallelogram
///
///```ignore
/// ____________ (P1)
/// / /
/// / . /
/// / /
/// /___________/ (P0)
/// ```
///
/// 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.
#[cfg_attr(feature = "serde", derive(Deserialize, Serialize))]
#[cfg_attr(feature = "fmt", derive(Debug))]
#[derive(Copy, Clone, Default, PartialEq)]
pub struct MBPgram<R> {
/// X coordinate at which the upward slope starts
/// X coordinate for point P0 in the diagram above.
pub b_start: R,
/// X coordinate at which the upward slope ends
/// X coordinate for point P1 in the diagram above.
pub b_end: R,
/// Vertical range of the graph
/// Y coordinate for point P1, negative Y coordinate for point P0, in the diagram above.
pub max_m: R,
}
/// evaluate `f(x)` at the provided `x`, where `f` is defined by these constraints:
/// 1. f(x0) = -ymax
/// 2. f(x1) = ymax
fn eval_bounded_line<R: Real>(x: R, x0: R, x1: R, ymax: R) -> R {
let x = x - x0;
let x1 = x1 - x0;
// Now, f(0) = -ymax, f(x1) = ymax
let x = x.max_or_undefined(R::zero()).min_or_undefined(x1);
ymax * (R::two()*x/x1 - R::one())
}
impl<R> MBPgram<R> {
pub fn new(b_start: R, b_end: R, max_m: R) -> Self {
Self { b_start, b_end, max_m }
}
}
impl<R: Real> MBPgram<R> {
/// Return the new `M`
pub fn move_b(self, m: R, target_b: R) -> R {
// 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_or_undefined(y_right).min_or_undefined(y_left)
}
}
impl<R: Real> Material<R> for MBPgram<R> {
fn move_b_vec(&self, m: Vec3<R>, target_b: Vec3<R>) -> Vec3<R> {
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()),
)
}
}
impl<R: Real, P: Peano> DiscriminantCodable<P> for MBPgram<R> {
fn decode_discr(&self) -> Discr<P> {
let max_m = self.max_m;
let d = if max_m < R::zero() {
(-max_m.to_f32()) as u32
} else {
0
};
Discr::new(d)
}
fn encode_discr(d: Discr<P>) -> Self {
Self {
max_m: R::from_primitive(-(d.value() as i32)),
..Default::default()
}
}
}
#[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);
}
}