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

Have the GenericMaterial be an enum to make it easier to replace the ferromagnet model

Browse Source
This commit is contained in:
Colin
2020-08-27 12:51:12 -07:00
parent c68964403d
commit 6ac0dafa47
4 changed files with 152 additions and 95 deletions
Download Patch File Download Diff File Expand all files Collapse all files

11
examples/em_reflection.rs
View File

@@ -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();
}
}

12
examples/ferromagnet.rs
View File

@@ -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();
}
}

209
src/lib.rs
View File

@@ -51,7 +51,7 @@ pub mod consts {
#[derive(Default)]
pub struct SimState {
cells: Array2<Cell<GenericMaterial>>,
cells: Array2<Cell<mat::GenericMaterial>>,
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<GenericMaterial> {
self.cells[[y, x]].clone()
pub fn get(&self, x: usize, y: usize) -> &Cell<mat::GenericMaterial> {
&self.cells[[y, x]]
}
pub fn get_mut(&mut self, x: usize, y: usize) -> &mut Cell<GenericMaterial> {
pub fn get_mut(&mut self, x: usize, y: usize) -> &mut Cell<mat::GenericMaterial> {
&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<M> {
state: CellState,
mat: M,
@@ -167,7 +167,7 @@ impl<M> Cell<M> {
}
}
impl<M: Material + Clone> Cell<M> {
impl<M: Material> Cell<M> {
pub fn bz(&self) -> f64 {
consts::MU0 * (self.hz() + self.mat.mz())
}
@@ -176,7 +176,7 @@ impl<M: Material + Clone> Cell<M> {
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<M: Material + Clone> Cell<M> {
/// 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<Conductor> for GenericMaterial {
fn from(f: Conductor) -> Self {
Self::Conductor(f)
}
}
impl From<ComsolFerromagnet> 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::*;

15
src/render.rs
View File

@@ -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");