fdtd-coremem/crates/coremem/src/mat/linear.rs

use crate::CellState;
use crate::geom::Vec3;
use crate::mat::Material;
use crate::real::Real;
use crate::sim::StepParametersMut;

use serde::{Serialize, Deserialize};

/// Material which has a conductivity parameter, but cannot be magnetized
#[derive(Copy, Clone, Default, PartialEq, Serialize, Deserialize)]
pub struct Conductor<V> {
    pub conductivity: V,
}

pub type IsomorphicConductor<R> = Conductor<R>;
pub type AnisomorphicConductor<R> = Conductor<Vec3<R>>;

impl<V> Conductor<V> {
    pub fn new(conductivity: V) -> Self {
        Self {
            conductivity
        }
    }
}

impl<R: Real> Into<AnisomorphicConductor<R>> for IsomorphicConductor<R> {
    fn into(self) -> AnisomorphicConductor<R> {
        AnisomorphicConductor::new(Vec3::uniform(self.conductivity))
    }
}

impl<R: Real> AnisomorphicConductor<R> {
    pub fn new_anisotropic<R2: Real>(c: Vec3<R2>) -> Self {
        Self {
            conductivity: Vec3::new(c.x().cast(), c.y().cast(), c.z().cast())
        }
    }
}

impl<R: Real> Material<R> for AnisomorphicConductor<R> {
    fn step_parameters_mut<'a>(&'a mut self) -> StepParametersMut<'a, R> {
        StepParametersMut::default().with_conductivity(self.conductivity)
    }
}

impl<R: Real> Material<R> for IsomorphicConductor<R> {
    fn step_parameters_mut<'a>(&'a mut self) -> StepParametersMut<'a, R> {
        StepParametersMut::default().with_conductivity(Vec3::uniform(self.conductivity))
    }
}

/// Material which can be magnetized, but has no hysteresis and no coercivity.
#[derive(Copy, Clone, Default, PartialEq, Serialize, Deserialize)]
pub struct LinearMagnet<R> {
    /// \mu_r
    relative_permeability: Vec3<R>,
    m: Vec3<R>,
}

impl<R: Real> LinearMagnet<R> {
    pub fn new<R2: Real>(relative_permeability: R2) -> Self {
        Self {
            relative_permeability: Vec3::uniform(relative_permeability).cast(),
            m: Vec3::zero(),
        }
    }
    pub fn new_anisotropic<R2: Real>(relative_permeability: Vec3<R2>) -> Self {
        Self {
            relative_permeability: relative_permeability.cast(),
            m: Vec3::zero()
        }
    }
}

impl<R: Real> Material<R> for LinearMagnet<R> {
    fn m(&self) -> Vec3<R> {
        self.m
    }
    fn step_b(&mut self, _context: &CellState<R>, delta_b: Vec3<R>) {
        //```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;
    }
}

#[cfg(test)]
mod test {
    use super::*;
    use float_eq::assert_float_eq;

    #[test]
    fn linear_magnet_steep() {
        let mut mag = LinearMagnet::<f64>::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::<f64>::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::<f32>::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);
    }
}