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

Model the H/B field in 3d, mostly

Browse Source
This commit is contained in:
Colin
2020-09-17 22:53:12 -07:00
parent 54b23a0494
commit a1f5b96ae1
5 changed files with 68 additions and 92 deletions
Download Patch File Download Diff File Expand all files Collapse all files

102
src/mat.rs
View File

@@ -1,5 +1,7 @@
use crate::{CellState, consts, flt::Flt};
use crate::{CellState, consts};
use crate::flt::Flt;
use crate::geom::{Line, Point, Polygon};
use crate::vec3::Vec3;
use enum_dispatch::enum_dispatch;
use std::cmp::Ordering;
@@ -10,21 +12,21 @@ pub trait Material {
fn conductivity(&self) -> Flt {
0.0
}
/// Return the magnetization in the z direction (M_z).
fn mz(&self) -> Flt {
0.0
/// Return the magnetization.
fn m(&self) -> Vec3 {
Vec3::zero()
}
/// Return the new h_z, and optionally change any internal state (e.g. magnetization).
fn step_h(&mut self, context: &CellState, delta_bz: Flt) -> Flt {
/// Return the new h, and optionally change any internal state (e.g. magnetization).
fn step_h(&mut self, context: &CellState, delta_b: Vec3) -> Vec3 {
// B = mu0*(H+M), and in a default material M=0.
context.hz() + delta_bz / consts::MU0
context.h() + delta_b / consts::real::MU0()
}
}
#[derive(Clone, Default)]
pub struct Static {
pub conductivity: Flt,
pub mz: Flt,
pub m: Vec3,
}
impl Static {
@@ -40,57 +42,8 @@ impl Material for Static {
fn conductivity(&self) -> Flt {
self.conductivity
}
fn mz(&self) -> Flt {
self.mz
}
}
#[derive(Clone, Default)]
pub struct ComsolFerromagnet {
/// Instantaneous magnetization in the z direction.
pub mz: Flt,
// Ferromagnetic parameters:
/// Magnetic saturation (symmetric, so saturates to either +Ms or -Ms).
pub ms: Flt,
/// Rate at which M changes w.r.t. H.
pub xi: Flt,
}
impl Material for ComsolFerromagnet {
fn mz(&self) -> Flt {
self.mz
}
fn step_h(&mut self, context: &CellState, delta_bz: Flt) -> Flt {
// COMSOL ferromagnetic model: https://www.comsol.com/blogs/accessing-external-material-models-for-magnetic-simulations/
// delta(M) = xi * delta*(H) # delta(M) = M(t+1) - M(t)
// new_M = clamp(M+delta(M), -Ms, Ms)) # i.e. clamp to saturation
// Where delta*(H) is H(t+1) - H(t), but clamping H(t) to <= 0 if positively saturated
// and clamping H(t) to >= 0 if negatively saturated
let expected_delta_h = delta_bz / consts::MU0;
let expected_new_h = context.hz() + expected_delta_h;
// let unsaturated_new_m = self.mz + self.xi * unsaturated_delta_h;
let unsaturated_new_m = if self.mz == self.ms {
// positively saturated
let delta_h = expected_new_h.min(0.0) - context.hz().min(0.0);
// assert!(self.ms == 0.0, "pos sat");
self.mz + self.xi * delta_h
} else if self.mz == -self.ms {
// negatively saturated
let delta_h = expected_new_h.max(0.0) - context.hz().max(0.0);
//assert!(self.ms == 0.0, "neg sat");
self.mz + self.xi * delta_h
} else {
// not saturated
self.mz + self.xi * expected_delta_h
};
let new_m = unsaturated_new_m.min(self.ms).max(-self.ms);
let delta_m = new_m - self.mz;
self.mz = new_m;
// B = mu0*(H + M)
// \Delta B = mu0*(\Delta H + \Delta M)
let delta_h = delta_bz / consts::MU0 - delta_m;
context.hz() + delta_h
fn m(&self) -> Vec3 {
self.m
}
}
@@ -169,9 +122,11 @@ impl MHCurve {
#[derive(Clone)]
pub struct PiecewiseLinearFerromagnet {
/// Instantaneous magnetization in the z direction.
pub mz: Flt,
pub m: Vec3,
pub conductivity: Flt,
mh_curve: MHCurve,
mh_curve_x: MHCurve,
mh_curve_y: MHCurve,
mh_curve_z: MHCurve,
}
impl PiecewiseLinearFerromagnet {
@@ -184,28 +139,33 @@ impl PiecewiseLinearFerromagnet {
}).collect();
Self {
mz: 0.0,
m: Vec3::default(),
conductivity: 0.0,
mh_curve: MHCurve::new(&*mh_points),
mh_curve_x: MHCurve::new(&*mh_points),
mh_curve_y: MHCurve::new(&*mh_points),
mh_curve_z: MHCurve::new(&*mh_points),
}
}
}
impl Material for PiecewiseLinearFerromagnet {
fn mz(&self) -> Flt {
self.mz
fn m(&self) -> Vec3 {
self.m
}
fn conductivity(&self) -> Flt {
self.conductivity
}
fn step_h(&mut self, context: &CellState, delta_bz: Flt) -> Flt {
let (h, m) = (context.hz(), self.mz);
let target_hm = h + m + delta_bz/consts::MU0;
fn step_h(&mut self, context: &CellState, delta_b: Vec3) -> Vec3 {
let (h, m) = (context.h(), self.m);
let target_hm = h + m + delta_b/consts::real::MU0();
let (h, m) = self.mh_curve.move_to(h, m, target_hm);
// TODO: this is probably not the best way to generalize a BH curve into 3d.
let (hx, mx) = self.mh_curve_x.move_to(h.x(), m.x(), target_hm.x());
let (hy, my) = self.mh_curve_x.move_to(h.y(), m.y(), target_hm.y());
let (hz, mz) = self.mh_curve_x.move_to(h.z(), m.z(), target_hm.z());
self.mz = m;
h
self.m = Vec3::new(mx, my, mz);
Vec3::new(hx, hy, hz)
}
}

17
src/meas.rs
View File

@@ -49,8 +49,8 @@ pub struct Magnetization(pub u32, pub u32);
impl AbstractMeasurement for Magnetization {
fn eval(&self, state: &SimSnapshot) -> String {
let mz = state.get(self.0, self.1).mat().mz();
format!("Mz({}, {}): {:.2e}", self.0, self.1, mz)
let m = state.get(self.0, self.1).mat().m();
format!("M({}, {}): ({:.2e}, {:.2e}, {:.2e})", self.0, self.1, m.x(), m.y(), m.z())
}
}
@@ -59,8 +59,8 @@ pub struct MagneticFlux(pub u32, pub u32);
impl AbstractMeasurement for MagneticFlux {
fn eval(&self, state: &SimSnapshot) -> String {
let bz = state.get(self.0, self.1).bz();
format!("Bz({}, {}): {:.2e}", self.0, self.1, bz)
let b = state.get(self.0, self.1).b();
format!("Bz({}, {}): ({:.2e}, {:.2e}, {:.2e})", self.0, self.1, b.x(), b.y(), b.z())
}
}
@@ -69,8 +69,8 @@ pub struct MagneticStrength(pub u32, pub u32);
impl AbstractMeasurement for MagneticStrength {
fn eval(&self, state: &SimSnapshot) -> String {
let bz = state.get(self.0, self.1).hz();
format!("Hz({}, {}): {:.2e}", self.0, self.1, bz)
let h = state.get(self.0, self.1).h();
format!("Hz({}, {}): ({:.2e}, {:.2e}, {:.2e})", self.0, self.1, h.x(), h.y(), h.z())
}
}
@@ -82,9 +82,10 @@ impl AbstractMeasurement for Energy {
#[allow(non_snake_case)]
let dV = f*f*f;
let e = state.map_sum(|cell| {
0.5 * cell.hz() * cell.bz() * dV
0.5 * cell.h().z() * cell.b().z() * dV
});
format!("U: {:.2e}", e)
// TODO: update to 3 dimensions!
format!("U(TODO): {:.2e}", e)
}
}

4
src/render.rs
View File

@@ -81,9 +81,9 @@ impl<'a> RenderSteps<'a> {
for y in 0..h {
for x in 0..w {
let cell = self.sim.get(x, y);
let r = scale_signed_to_u8(cell.mat().mz(), 100.0);
let r = scale_signed_to_u8(cell.mat().m().z(), 100.0);
let b = scale_unsigned_to_u8(cell.mat().conductivity(), 10.0);
let g = scale_signed_to_u8(cell.bz(), 1.0e-4);
let g = scale_signed_to_u8(cell.b().z(), 1.0e-4);
self.im.put_pixel(x, y, Rgb([r, g, b]));
}
}

31
src/sim.rs
View File

@@ -75,14 +75,14 @@ impl<M: Material + Default + Send + Sync> SimState<M> {
let mut cell: Cell<mat::Static> = Default::default();
for sample in cells {
cell.state.e += sample.state.e();
cell.state.hz += sample.hz();
cell.state.h += sample.h();
cell.mat.conductivity += sample.mat.conductivity();
cell.mat.mz += sample.mat.mz();
cell.mat.m += sample.mat.m();
}
cell.state.e /= Real::from_inner(sample_area);
cell.state.hz /= sample_area;
cell.state.h /= Real::from_inner(sample_area);
cell.mat.conductivity /= sample_area;
cell.mat.mz /= sample_area;
cell.mat.m /= Real::from_inner(sample_area);
cell
});
@@ -181,6 +181,9 @@ impl<M> Cell<M> {
pub fn hz(&self) -> Flt {
self.state.hz()
}
pub fn h(&self) -> Vec3 {
self.state.h()
}
pub fn mat(&self) -> &M {
&self.mat
}
@@ -190,8 +193,11 @@ impl<M> Cell<M> {
}
impl<M: Material> Cell<M> {
pub fn b(&self) -> Vec3 {
(self.h() + self.mat.m()) * consts::real::MU0()
}
pub fn bz(&self) -> Flt {
consts::MU0 * (self.hz() + self.mat.mz())
self.b().z()
}
pub fn current(&self) -> Point {
@@ -200,7 +206,7 @@ impl<M: Material> Cell<M> {
}
fn impulse_bz(&mut self, delta_bz: Flt) {
self.state.hz = self.mat.step_h(&self.state, delta_bz).into();
self.state.h = self.mat.step_h(&self.state, Vec3::new(0.0, 0.0, delta_bz));
}
fn step_h(mut self, right: &Self, down: &Self, delta_t: Real, feature_size: Real) -> Self {
@@ -227,10 +233,10 @@ impl<M: Material> Cell<M> {
let nabla_e = Vec3::from((delta_ez_y - delta_ey_z, delta_ex_z - delta_ez_x, delta_ey_x - delta_ex_y));
let delta_b = -nabla_e * delta_t;
let hz = self.mat.step_h(&self.state, delta_b.z());
let h = self.mat.step_h(&self.state, delta_b);
Cell {
state: CellState {
hz: hz.into(),
h: h.into(),
..self.state
},
mat: self.mat,
@@ -284,7 +290,7 @@ impl<M: Material> Cell<M> {
Cell {
state: CellState {
e: e_next,
hz: self.state.hz,
h: self.state.h,
},
mat: self.mat,
}
@@ -294,7 +300,7 @@ impl<M: Material> Cell<M> {
#[derive(Copy, Clone, Default)]
pub struct CellState {
e: Vec3,
hz: Real,
h: Vec3,
}
impl CellState {
@@ -305,7 +311,10 @@ impl CellState {
self.e.y()
}
pub fn hz(&self) -> Flt {
self.hz.into()
self.h.z().into()
}
pub fn h(&self) -> Vec3 {
self.h
}
pub fn e(&self) -> Vec3 {
self.e

6
src/vec3.rs
View File

@@ -17,6 +17,12 @@ impl Vec3 {
z: z.into()
}
}
pub fn unit() -> Self {
Self::new(1.0, 1.0, 1.0)
}
pub fn zero() -> Self {
Self::new(0.0, 0.0, 0.0)
}
pub fn x(&self) -> Flt {
self.x.into()
}