diff --git a/examples/em_reflection.rs b/examples/em_reflection.rs index b082e26..5c5fdcd 100644 --- a/examples/em_reflection.rs +++ b/examples/em_reflection.rs @@ -1,6 +1,5 @@ -use coremem::SimState; +use coremem::{consts, mat, SimState}; use coremem::render::ColorTermRenderer as Renderer; -use coremem::consts; use std::{thread, time}; fn main() { @@ -9,11 +8,11 @@ fn main() { for inset in 0..20 { for x in 0..width { - state.get_mut(x, inset).mat_mut().conductivity += (0.1*(20.0 - inset as f64)); + *state.get_mut(x, inset).mat_mut() = mat::Conductor { conductivity: 0.1*(20.0 - inset as f64) }.into(); } for y in 0..100 { - state.get_mut(inset, y).mat_mut().conductivity += (0.1*(20.0 - inset as f64)); - state.get_mut(width-1 - inset, y).mat_mut().conductivity += (0.1*(20.0 - inset as f64)); + *state.get_mut(inset, y).mat_mut() = mat::Conductor { conductivity: 0.1*(20.0 - inset as f64) }.into(); + *state.get_mut(width-1 - inset, y).mat_mut() = mat::Conductor { conductivity: 0.1*(20.0 - inset as f64) }.into(); } } for y in 75..100 { @@ -22,7 +21,7 @@ fn main() { // NB: different sources give pretty different values for this // NB: Simulation misbehaves for values > 10... Proably this model isn't so great. // Maybe use \eps or \xi instead of conductivity. - state.get_mut(x, y).mat_mut().conductivity = (2.17).into(); + *state.get_mut(x, y).mat_mut() = mat::Conductor { conductivity: 2.17 }.into(); } } diff --git a/examples/ferromagnet.rs b/examples/ferromagnet.rs index df8b48a..75ffb8b 100644 --- a/examples/ferromagnet.rs +++ b/examples/ferromagnet.rs @@ -1,6 +1,5 @@ -use coremem::SimState; +use coremem::{consts, mat, SimState}; use coremem::render::ColorTermRenderer as Renderer; -use coremem::consts; use std::{thread, time}; fn main() { @@ -9,13 +8,16 @@ fn main() { for y in 0..100 { for x in 50..60 { - state.get_mut(x, y).mat_mut().conductivity = 10.0; + *state.get_mut(x, y).mat_mut() = mat::Conductor { conductivity: 10.0 }.into(); } } for y in 40..60 { for x in 62..70 { - state.get_mut(x, y).mat_mut().xi = 150.0; - state.get_mut(x, y).mat_mut().ms = 2.0; + *state.get_mut(x, y).mat_mut() = mat::ComsolFerromagnet { + xi: 150.0, + ms: 2.0, + ..mat::ComsolFerromagnet::default() + }.into(); } } diff --git a/src/lib.rs b/src/lib.rs index 7834a3f..3b2d7c3 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -51,7 +51,7 @@ pub mod consts { #[derive(Default)] pub struct SimState { - cells: Array2>, + cells: Array2>, feature_size: R64, step_no: u64, } @@ -79,7 +79,7 @@ impl SimState { let cell = self.get(right_x-1, down_y-1); let right_cell = self.get(right_x, down_y-1); let down_cell = self.get(right_x-1, down_y); - working_cells[[down_y-1, right_x-1]] = cell.step_h(right_cell, down_cell, half_time_step.into()); + working_cells[[down_y-1, right_x-1]] = cell.clone().step_h(right_cell, down_cell, half_time_step.into()); } } std::mem::swap(&mut working_cells, &mut self.cells); @@ -90,7 +90,7 @@ impl SimState { let cell = self.get(x, y); let left_cell = self.get(x-1, y); let up_cell = self.get(x, y-1); - working_cells[[y, x]] = cell.step_e(left_cell, up_cell, half_time_step.into(), self.feature_size.into()); + working_cells[[y, x]] = cell.clone().step_e(left_cell, up_cell, half_time_step.into(), self.feature_size.into()); } } std::mem::swap(&mut working_cells, &mut self.cells); @@ -113,10 +113,10 @@ impl SimState { pub fn height(&self) -> usize { self.cells.shape()[0] } - pub fn get(&self, x: usize, y: usize) -> Cell { - self.cells[[y, x]].clone() + pub fn get(&self, x: usize, y: usize) -> &Cell { + &self.cells[[y, x]] } - pub fn get_mut(&mut self, x: usize, y: usize) -> &mut Cell { + pub fn get_mut(&mut self, x: usize, y: usize) -> &mut Cell { &mut self.cells[[y, x]] } @@ -143,7 +143,7 @@ impl SimState { /// The dot on bottom is Ex of the cell at (x, y+1) and the dot on the right is the Ey of the cell at /// (x+1, y). The `+` only indicates the corner of the cell -- nothing of interest is measured at /// the pluses. -#[derive(Copy, Clone, Default)] +#[derive(Clone, Default)] pub struct Cell { state: CellState, mat: M, @@ -167,7 +167,7 @@ impl Cell { } } -impl Cell { +impl Cell { pub fn bz(&self) -> f64 { consts::MU0 * (self.hz() + self.mat.mz()) } @@ -176,7 +176,7 @@ impl Cell { self.state.hz = self.mat.step_h(&self.state, delta_bz).into(); } - fn step_h(mut self, right: Self, down: Self, _delta_t: R64) -> Self { + fn step_h(mut self, right: &Self, down: &Self, _delta_t: R64) -> Self { // ```tex // Maxwell's equation: $\nabla x E = -dB/dt$ // Expand: $dE_y/dx - dE_x/dy = -dB_z/dt$ @@ -202,7 +202,7 @@ impl Cell { /// delta_t = timestep covered by this function. i.e. it should be half the timestep of the simulation /// feature_size = how many units apart is the center of each adjacent cell on the grid. - fn step_e(self, left: Self, up: Self, delta_t: R64, feature_size: R64) -> Self { + fn step_e(self, left: &Self, up: &Self, delta_t: R64, feature_size: R64) -> Self { // ```tex // Ampere's circuital law with Maxwell's addition, in SI units: // $\nabla x B = \mu_0 (J + \epsilon_0 dE/dt)$ where J = current density = $\sigma E$, $\sigma$ being a material parameter @@ -268,74 +268,131 @@ impl CellState { } } -pub trait Material { - /// Return \sigma, the electrical conductivity. - /// For a vacuum, this is zero. For a perfect conductor, \inf. - fn conductivity(&self) -> f64 { - 0.0 +pub mod mat { + use super::{CellState, consts}; + pub trait Material { + /// Return \sigma, the electrical conductivity. + /// For a vacuum, this is zero. For a perfect conductor, \inf. + fn conductivity(&self) -> f64 { + 0.0 + } + /// Return the magnetization in the z direction (M_z). + fn mz(&self) -> f64 { + 0.0 + } + /// Return the new h_z, and optionally change any internal state (e.g. magnetization). + fn step_h(&mut self, context: &CellState, delta_bz: f64) -> f64 { + // B = mu0*(H+M), and in a default material M=0. + context.hz() + delta_bz / consts::MU0 + } } - /// Return the magnetization in the z direction (M_z). - fn mz(&self) -> f64 { - 0.0 + + #[derive(Clone, Default)] + pub struct Conductor { + pub conductivity: f64, } - /// Return the new h_z, and optionally change any internal state (e.g. magnetization). - fn step_h(&mut self, context: &CellState, delta_bz: f64) -> f64 { - // B = mu0*(H+M), and in a default material M=0. - context.hz() + delta_bz / consts::MU0 + + impl Material for Conductor { + fn conductivity(&self) -> f64 { + self.conductivity + } + } + + #[derive(Clone, Default)] + pub struct ComsolFerromagnet { + /// Instantaneous magnetization in the z direction. + pub mz: f64, + // Ferromagnetic parameters: + /// Magnetic saturation (symmetric, so saturates to either +Ms or -Ms). + pub ms: f64, + /// Rate at which M changes w.r.t. H. + pub xi: f64, + } + + impl Material for ComsolFerromagnet { + fn mz(&self) -> f64 { + self.mz + } + fn step_h(&mut self, context: &CellState, delta_bz: f64) -> f64 { + // 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 + } + } + + #[derive(Clone)] + pub enum GenericMaterial { + Conductor(Conductor), + ComsolFerromagnet(ComsolFerromagnet), + } + + impl From for GenericMaterial { + fn from(f: Conductor) -> Self { + Self::Conductor(f) + } + } + + impl From for GenericMaterial { + fn from(f: ComsolFerromagnet) -> Self { + Self::ComsolFerromagnet(f) + } + } + + impl Default for GenericMaterial { + fn default() -> Self { + Conductor::default().into() + } + } + + impl Material for GenericMaterial { + fn conductivity(&self) -> f64 { + match self { + Self::Conductor(inner) => inner.conductivity(), + Self::ComsolFerromagnet(inner) => inner.conductivity(), + } + } + /// Return the magnetization in the z direction (M_z). + fn mz(&self) -> f64 { + match self { + Self::Conductor(inner) => inner.mz(), + Self::ComsolFerromagnet(inner) => inner.mz(), + } + } + /// Return the new h_z, and optionally change any internal state (e.g. magnetization). + fn step_h(&mut self, context: &CellState, delta_bz: f64) -> f64 { + match self { + Self::Conductor(inner) => inner.step_h(context, delta_bz), + Self::ComsolFerromagnet(inner) => inner.step_h(context, delta_bz), + } + } } } -#[derive(Clone, Default)] -pub struct GenericMaterial { - // Electrical parameters: - pub conductivity: f64, - // Magnetic state: - /// Instantaneous magnetization in the z direction. - pub mz: f64, - // Ferromagnetic parameters: - /// Magnetic saturation (symmetric, so saturates to either +Ms or -Ms). - pub ms: f64, - /// Rate at which M changes w.r.t. H. - pub xi: f64, -} - -impl Material for GenericMaterial { - fn conductivity(&self) -> f64 { - self.conductivity - } - fn mz(&self) -> f64 { - self.mz - } - fn step_h(&mut self, context: &CellState, delta_bz: f64) -> f64 { - // 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 - } -} +pub use mat::*; diff --git a/src/render.rs b/src/render.rs index f6c1d75..487db32 100644 --- a/src/render.rs +++ b/src/render.rs @@ -1,6 +1,5 @@ use ansi_term::Color::RGB; -use crate::consts; -use crate::SimState; +use crate::{consts, Material as _, SimState}; use std::fmt::Write as _; pub struct NumericTermRenderer; @@ -33,9 +32,9 @@ fn norm_color(v: f64) -> u8 { (v * 64.0 + 128.0).max(0.0).min(255.0) as u8 } -fn curl(x: f32, y: f32) -> f32 { +fn curl(x: f64, y: f64) -> f64 { let c = x * y; - if c >= 0f32 { + if c >= 0.0 { c.sqrt() } else { -(-c).sqrt() @@ -51,14 +50,14 @@ impl ColorTermRenderer { let cell = state.get(x, y); //let r = norm_color(cell.bz() * consts::C); //let r = 0; - let r = norm_color(cell.mat.mz); - let b = (55.0*cell.mat().conductivity).min(255.0) as u8; + let r = norm_color(cell.mat.mz()); + let b = (55.0*cell.mat().conductivity()).min(255.0) as u8; //let b = 0; //let g = norm_color(cell.ex()); //let b = norm_color(cell.ey()); //let g = norm_color(curl(cell.ex(), cell.ey())); - //let g = norm_color((cell.bz() * 3.0e8).into()); - let g = norm_color(cell.ey().into()); + let g = norm_color((cell.bz() * 3.0e8).into()); + //let g = norm_color(cell.ey().into()); write!(&mut buf, "{}", RGB(r, g, b).paint(square)); } write!(&mut buf, "\n");