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

Factor mat.rs and sim.rs out of lib.rs

Browse Source
This commit is contained in:
Colin
2020-08-27 19:17:07 -07:00
parent 57258eb6cc
commit b25ccc2471
4 changed files with 395 additions and 389 deletions
Download Patch File Download Diff File Expand all files Collapse all files

393
src/lib.rs
View File

@@ -6,9 +6,13 @@
//! [1] https://www.eecs.wsu.edu/~schneidj/ufdtd/ufdtd.pdf
use decorum::R64;
use ndarray::Array2;
pub mod mat;
pub mod render;
pub mod sim;
pub use mat::*;
pub use sim::*;
// Some things to keep in mind:
// B = mu_r*H + M
@@ -52,390 +56,3 @@ pub mod consts {
}
}
#[derive(Default)]
pub struct SimState {
cells: Array2<Cell<mat::GenericMaterial>>,
feature_size: R64,
step_no: u64,
}
impl SimState {
pub fn new(width: usize, height: usize, feature_size: f64) -> Self {
Self {
cells: Array2::default((height, width)),
feature_size: feature_size.into(),
..Default::default()
}
}
pub fn time(&self) -> f64 {
(self.timestep() * self.step_no as f64).into()
}
pub fn step(&mut self) {
use consts::real::*;
let half_time_step = HALF() * self.timestep();
let mut working_cells = Array2::default((self.height(), self.width()));
// first advance all the magnetic fields
for down_y in 1..self.height() {
for right_x in 1..self.width() {
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.clone().step_h(right_cell, down_cell, half_time_step.into());
}
}
std::mem::swap(&mut working_cells, &mut self.cells);
// now advance electic fields
for y in 1..self.height() {
for x in 1..self.width() {
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.clone().step_e(left_cell, up_cell, half_time_step.into(), self.feature_size.into());
}
}
std::mem::swap(&mut working_cells, &mut self.cells);
self.step_no += 1;
}
pub fn impulse_ex(&mut self, x: usize, y: usize, ex: f64) {
self.cells[[y, x]].state.ex += ex;
}
pub fn impulse_ey(&mut self, x: usize, y: usize, ey: f64) {
self.cells[[y, x]].state.ey += ey;
}
pub fn impulse_bz(&mut self, x: usize, y: usize, bz: f64) {
self.cells[[y, x]].impulse_bz(bz);
}
pub fn width(&self) -> usize {
self.cells.shape()[1]
}
pub fn height(&self) -> usize {
self.cells.shape()[0]
}
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<mat::GenericMaterial> {
&mut self.cells[[y, x]]
}
fn timestep(&self) -> R64 {
self.feature_size / consts::real::C()
}
}
/// Conceptually, one cell looks like this:
///
/// +-------.-------+
/// | Ex |
/// | |
/// | |
/// .Ey .Bz .
/// | |
/// | |
/// | |
/// +-------.-------+
///
/// Where the right hand rule indicates that positive Bz is pointing into the page, away from the
/// reader.
///
/// 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(Clone, Default)]
pub struct Cell<M> {
state: CellState,
mat: M,
}
impl<M> Cell<M> {
pub fn ex(&self) -> f64 {
self.state.ex()
}
pub fn ey(&self) -> f64 {
self.state.ey()
}
pub fn hz(&self) -> f64 {
self.state.hz()
}
pub fn mat(&self) -> &M {
&self.mat
}
pub fn mat_mut(&mut self) -> &mut M {
&mut self.mat
}
}
impl<M: Material> Cell<M> {
pub fn bz(&self) -> f64 {
consts::MU0 * (self.hz() + self.mat.mz())
}
fn bz_to_hz(&self, bz: f64) -> f64 {
// B = mu0*(H + M) => H = B/mu0 - M
bz/consts::MU0 - self.mat.mz()
}
fn impulse_bz(&mut self, delta_bz: f64) {
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 {
// ```tex
// Maxwell-Faraday equation: $\nabla x E = -dB/dt$
// Expand: $dE_y/dx - dE_x/dy = -dB_z/dt$
// Rearrange: $dB_z/dt = dE_x/dy - dE_y/dx$
// Discretize: $(\Delta B_z)/(\Delta t) = (\Delta E_x)/(\Delta y) - (\Delta E_y)/(\Delta x)$
//
// light travels C meters per second, and we're only stepping half a timestep here (the other half when we advance E),
// therefore $(\Delta x)/(2 \Delta t) = (\Delta y)/(2 \Delta t) = c$ if we use SI units.
// Simplify: $\Delta B_z = (\Delta E_x)/(2c) - (\Delta E_y)/(2c)$
// ```
let delta_ex = down.ex() - self.ex();
let delta_ey = right.ey() - self.ey();
let delta_bz = (delta_ex - delta_ey) / (2.0*consts::C);
let hz = self.mat.step_h(&self.state, delta_bz);
Cell {
state: CellState {
hz: hz.into(),
..self.state
},
mat: self.mat,
}
}
/// 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 {
// ```tex
// Ampere's circuital law with Maxwell's addition, in SI units ("macroscopic version"):
// $\nabla x H = J_f + dD/dt$ where $J_f$ = current density = $\sigma E$, $\sigma$ being a material parameter ("conductivity")
// Note that $D = \epsilon_0 E + P$, but we don't simulate any material where $P \ne 0$.
// Expand $D$: $\nabla x H = J_f + \epsilon_0 dE/dt$
// Expand $J$: $\nabla x H = \sigma E + \epsilon_0 dE/dt$
// Rearrange: $dE/dt = 1/\epsilon_0 (\nabla x H - \sigma E)$
// Expand: $dE_x/dt = 1/\epsilon_0 (dH_z/dy - \sigma E_x)$ (1); $dE_y/dt = 1/\epsilon_0 (-dH_z/dx - \sigma E_y)$ (2)
//
// Consider (1): let $E_p$ be $E_x$ at $T-\Delta t$ and $E_n$ be $E_x$ at $T+\Delta t$.
// Linear expansion about $t=T$, and discretized:
// $(E_n-E_p)/(2\Delta{t}) = 1/\epsilon_0(\Delta{H_z}/\Delta{y} - \sigma(E_n+E_p)/2)$
// Normalize: $E_n - E_p = 2\Delta{t}/\epsilon_0 \Delta{H_z}/\Delta{y} - \sigma/\epsilon_0 \Delta{t} (E_n + E_p)$
// Rearrange: $E_n(1 + \sigma/\epsilon_0 \Delta{t}) = E_p(1 - \sigma/\epsilon_0 \Delta{t}) + 2\Delta{t}/\epsilon_0 \Delta{H_z}/\Delta{y}$
// Then $E_n$ (i.e. the x value of $E$ after this step) is trivially solved
//
// Consider (2): let $E_p$ be $E_y$ at $T-\Delta t$ and $E_n$ be $E_y$ at $T+\Delta t$.
// Linear expansion about $t=T$, and discretized:
// $(E_n-E_p)/(2\Delta{t}) = 1/\epsilon_0(-\Delta{H_z}/\Delta{x} - \sigma(E_n+E_p)/2)$
// Normalize: $E_n - E_p = -2\Delta{t}/\epsilon_0 \Delta{H_z}/\Delta{x} - \sigma/\epsilon_0 \Delta{t} (E_n + E_p)$
// Rearrange: $E_n(1 + \sigma/\epsilon_0 \Delta{t}) = E_p(1 - \sigma/\epsilon_0 \Delta{t}) - 2\Delta{t}/\epsilon_0 \Delta{H_z}/\Delta{x}$
// Then $E_n$ (i.e. the y value of $E$ after this step) is trivially solved
// ```
use consts::real::{C, C2, EPS0, HALF, ONE, TWO};
let sigma: R64 = self.mat.conductivity().into();
let delta_hz_y = self.hz() - self.bz_to_hz(up.bz());
let ex_rhs = self.state.ex*(ONE() - sigma/EPS0()*delta_t) + TWO()*delta_t/EPS0()*delta_hz_y/feature_size;
let ex_next = ex_rhs / (ONE() + sigma/EPS0()*delta_t);
let delta_hz_x = self.hz() - self.bz_to_hz(left.bz());
let ey_rhs = self.state.ey*(ONE() - sigma/EPS0()*delta_t) - TWO()*delta_t/EPS0()*delta_hz_x/feature_size;
let ey_next = ey_rhs / (ONE() + sigma/EPS0()*delta_t);
Cell {
state: CellState {
ex: ex_next,
ey: ey_next,
hz: self.state.hz,
},
mat: self.mat,
}
}
}
#[derive(Copy, Clone, Default)]
pub struct CellState {
ex: R64,
ey: R64,
hz: R64,
}
impl CellState {
pub fn ex(&self) -> f64 {
self.ex.into()
}
pub fn ey(&self) -> f64 {
self.ey.into()
}
pub fn hz(&self) -> f64 {
self.hz.into()
}
}
pub mod mat {
use super::{CellState, consts};
use enum_dispatch::enum_dispatch;
use piecewise_linear::PiecewiseLinearFunction;
#[enum_dispatch]
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
}
}
#[derive(Clone, Default)]
pub struct Conductor {
pub conductivity: f64,
}
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
}
}
pub struct BHCurve {
/// describes the curve as B rises from 0.
rising: PiecewiseLinearFunction<f64>,
/// describes the curve as B falls towards 0.
falling: PiecewiseLinearFunction<f64>,
}
impl BHCurve {
fn b_bounds_for_h(&self, h: f64) -> (f64, f64) {
let (rise_start, rise_end) = self.rising.domain();
let bottom = match h {
x if (..-rise_end).contains(&x) => {
// negative saturation
-self.falling.y_at_x(-rise_end).unwrap() + consts::MU0 * (h + rise_end)
},
x if (-rise_end..rise_start).contains(&x) => {
// negative
-self.falling.y_at_x(-h).unwrap()
}
x if (rise_start..rise_end).contains(&x) => {
// positive
self.rising.y_at_x(h).unwrap()
},
x if (rise_end..).contains(&x) => {
// positive saturation
self.rising.y_at_x(rise_end).unwrap() + consts::MU0 * (h - rise_end)
}
_ => unreachable!("invalid h: {}", h),
};
let top = match h {
x if (..-rise_end).contains(&x) => {
// negative saturation
-self.rising.y_at_x(-rise_end).unwrap() + consts::MU0 * (h + rise_end)
},
x if (-rise_end..-rise_start).contains(&x) => {
// negative
-self.rising.y_at_x(-h).unwrap()
},
x if (-rise_start..rise_end).contains(&x) => {
// positive
self.falling.y_at_x(h).unwrap()
},
x if (rise_end..).contains(&x) => {
// positive saturation
self.falling.y_at_x(rise_end).unwrap() + consts::MU0 * (h - rise_end)
},
_ => unreachable!("invalid h: {}", h),
};
(bottom, top)
}
}
pub struct BHFerromagnet {
/// Instantaneous magnetization in the z direction.
pub mz: f64,
bh_curve: BHCurve,
}
impl Material for BHFerromagnet {
fn step_h(&mut self, context: &CellState, delta_bz: f64) -> f64 {
unimplemented!()
}
}
#[enum_dispatch(Material)]
#[derive(Clone)]
pub enum GenericMaterial {
Conductor(Conductor),
ComsolFerromagnet(ComsolFerromagnet),
}
impl Default for GenericMaterial {
fn default() -> Self {
Conductor::default().into()
}
}
}
pub use mat::*;

158
src/mat.rs Normal file
View File

@@ -0,0 +1,158 @@
use crate::{CellState, consts};
use enum_dispatch::enum_dispatch;
use piecewise_linear::PiecewiseLinearFunction;
#[enum_dispatch]
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
}
}
#[derive(Clone, Default)]
pub struct Conductor {
pub conductivity: f64,
}
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
}
}
pub struct BHCurve {
/// describes the curve as B rises from 0.
rising: PiecewiseLinearFunction<f64>,
/// describes the curve as B falls towards 0.
falling: PiecewiseLinearFunction<f64>,
}
impl BHCurve {
fn b_bounds_for_h(&self, h: f64) -> (f64, f64) {
let (rise_start, rise_end) = self.rising.domain();
let bottom = match h {
x if (..-rise_end).contains(&x) => {
// negative saturation
-self.falling.y_at_x(-rise_end).unwrap() + consts::MU0 * (h + rise_end)
},
x if (-rise_end..rise_start).contains(&x) => {
// negative
-self.falling.y_at_x(-h).unwrap()
}
x if (rise_start..rise_end).contains(&x) => {
// positive
self.rising.y_at_x(h).unwrap()
},
x if (rise_end..).contains(&x) => {
// positive saturation
self.rising.y_at_x(rise_end).unwrap() + consts::MU0 * (h - rise_end)
}
_ => unreachable!("invalid h: {}", h),
};
let top = match h {
x if (..-rise_end).contains(&x) => {
// negative saturation
-self.rising.y_at_x(-rise_end).unwrap() + consts::MU0 * (h + rise_end)
},
x if (-rise_end..-rise_start).contains(&x) => {
// negative
-self.rising.y_at_x(-h).unwrap()
},
x if (-rise_start..rise_end).contains(&x) => {
// positive
self.falling.y_at_x(h).unwrap()
},
x if (rise_end..).contains(&x) => {
// positive saturation
self.falling.y_at_x(rise_end).unwrap() + consts::MU0 * (h - rise_end)
},
_ => unreachable!("invalid h: {}", h),
};
(bottom, top)
}
}
pub struct BHFerromagnet {
/// Instantaneous magnetization in the z direction.
pub mz: f64,
bh_curve: BHCurve,
}
impl Material for BHFerromagnet {
fn step_h(&mut self, context: &CellState, delta_bz: f64) -> f64 {
unimplemented!()
}
}
#[enum_dispatch(Material)]
#[derive(Clone)]
pub enum GenericMaterial {
Conductor(Conductor),
ComsolFerromagnet(ComsolFerromagnet),
}
impl Default for GenericMaterial {
fn default() -> Self {
Conductor::default().into()
}
}

2
src/render.rs
View File

@@ -50,7 +50,7 @@ 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 r = norm_color(cell.mat().mz());
let b = (55.0*cell.mat().conductivity()).min(255.0) as u8;
//let b = 0;
//let b = norm_color(cell.ey());

231
src/sim.rs Normal file
View File

@@ -0,0 +1,231 @@
use crate::{consts};
use crate::mat::{GenericMaterial, Material};
use decorum::R64;
use ndarray::Array2;
#[derive(Default)]
pub struct SimState {
cells: Array2<Cell<GenericMaterial>>,
feature_size: R64,
step_no: u64,
}
impl SimState {
pub fn new(width: usize, height: usize, feature_size: f64) -> Self {
Self {
cells: Array2::default((height, width)),
feature_size: feature_size.into(),
..Default::default()
}
}
pub fn time(&self) -> f64 {
(self.timestep() * self.step_no as f64).into()
}
pub fn step(&mut self) {
use consts::real::*;
let half_time_step = HALF() * self.timestep();
let mut working_cells = Array2::default((self.height(), self.width()));
// first advance all the magnetic fields
for down_y in 1..self.height() {
for right_x in 1..self.width() {
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.clone().step_h(right_cell, down_cell, half_time_step.into());
}
}
std::mem::swap(&mut working_cells, &mut self.cells);
// now advance electic fields
for y in 1..self.height() {
for x in 1..self.width() {
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.clone().step_e(left_cell, up_cell, half_time_step.into(), self.feature_size.into());
}
}
std::mem::swap(&mut working_cells, &mut self.cells);
self.step_no += 1;
}
pub fn impulse_ex(&mut self, x: usize, y: usize, ex: f64) {
self.cells[[y, x]].state.ex += ex;
}
pub fn impulse_ey(&mut self, x: usize, y: usize, ey: f64) {
self.cells[[y, x]].state.ey += ey;
}
pub fn impulse_bz(&mut self, x: usize, y: usize, bz: f64) {
self.cells[[y, x]].impulse_bz(bz);
}
pub fn width(&self) -> usize {
self.cells.shape()[1]
}
pub fn height(&self) -> usize {
self.cells.shape()[0]
}
pub fn get(&self, x: usize, y: usize) -> &Cell<GenericMaterial> {
&self.cells[[y, x]]
}
pub fn get_mut(&mut self, x: usize, y: usize) -> &mut Cell<GenericMaterial> {
&mut self.cells[[y, x]]
}
fn timestep(&self) -> R64 {
self.feature_size / consts::real::C()
}
}
/// Conceptually, one cell looks like this:
///
/// +-------.-------+
/// | Ex |
/// | |
/// | |
/// .Ey .Bz .
/// | |
/// | |
/// | |
/// +-------.-------+
///
/// Where the right hand rule indicates that positive Bz is pointing into the page, away from the
/// reader.
///
/// 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(Clone, Default)]
pub struct Cell<M> {
state: CellState,
mat: M,
}
impl<M> Cell<M> {
pub fn ex(&self) -> f64 {
self.state.ex()
}
pub fn ey(&self) -> f64 {
self.state.ey()
}
pub fn hz(&self) -> f64 {
self.state.hz()
}
pub fn mat(&self) -> &M {
&self.mat
}
pub fn mat_mut(&mut self) -> &mut M {
&mut self.mat
}
}
impl<M: Material> Cell<M> {
pub fn bz(&self) -> f64 {
consts::MU0 * (self.hz() + self.mat.mz())
}
fn bz_to_hz(&self, bz: f64) -> f64 {
// B = mu0*(H + M) => H = B/mu0 - M
bz/consts::MU0 - self.mat.mz()
}
fn impulse_bz(&mut self, delta_bz: f64) {
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 {
// ```tex
// Maxwell-Faraday equation: $\nabla x E = -dB/dt$
// Expand: $dE_y/dx - dE_x/dy = -dB_z/dt$
// Rearrange: $dB_z/dt = dE_x/dy - dE_y/dx$
// Discretize: $(\Delta B_z)/(\Delta t) = (\Delta E_x)/(\Delta y) - (\Delta E_y)/(\Delta x)$
//
// light travels C meters per second, and we're only stepping half a timestep here (the other half when we advance E),
// therefore $(\Delta x)/(2 \Delta t) = (\Delta y)/(2 \Delta t) = c$ if we use SI units.
// Simplify: $\Delta B_z = (\Delta E_x)/(2c) - (\Delta E_y)/(2c)$
// ```
let delta_ex = down.ex() - self.ex();
let delta_ey = right.ey() - self.ey();
let delta_bz = (delta_ex - delta_ey) / (2.0*consts::C);
let hz = self.mat.step_h(&self.state, delta_bz);
Cell {
state: CellState {
hz: hz.into(),
..self.state
},
mat: self.mat,
}
}
/// 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 {
// ```tex
// Ampere's circuital law with Maxwell's addition, in SI units ("macroscopic version"):
// $\nabla x H = J_f + dD/dt$ where $J_f$ = current density = $\sigma E$, $\sigma$ being a material parameter ("conductivity")
// Note that $D = \epsilon_0 E + P$, but we don't simulate any material where $P \ne 0$.
// Expand $D$: $\nabla x H = J_f + \epsilon_0 dE/dt$
// Expand $J$: $\nabla x H = \sigma E + \epsilon_0 dE/dt$
// Rearrange: $dE/dt = 1/\epsilon_0 (\nabla x H - \sigma E)$
// Expand: $dE_x/dt = 1/\epsilon_0 (dH_z/dy - \sigma E_x)$ (1); $dE_y/dt = 1/\epsilon_0 (-dH_z/dx - \sigma E_y)$ (2)
//
// Consider (1): let $E_p$ be $E_x$ at $T-\Delta t$ and $E_n$ be $E_x$ at $T+\Delta t$.
// Linear expansion about $t=T$, and discretized:
// $(E_n-E_p)/(2\Delta{t}) = 1/\epsilon_0(\Delta{H_z}/\Delta{y} - \sigma(E_n+E_p)/2)$
// Normalize: $E_n - E_p = 2\Delta{t}/\epsilon_0 \Delta{H_z}/\Delta{y} - \sigma/\epsilon_0 \Delta{t} (E_n + E_p)$
// Rearrange: $E_n(1 + \sigma/\epsilon_0 \Delta{t}) = E_p(1 - \sigma/\epsilon_0 \Delta{t}) + 2\Delta{t}/\epsilon_0 \Delta{H_z}/\Delta{y}$
// Then $E_n$ (i.e. the x value of $E$ after this step) is trivially solved
//
// Consider (2): let $E_p$ be $E_y$ at $T-\Delta t$ and $E_n$ be $E_y$ at $T+\Delta t$.
// Linear expansion about $t=T$, and discretized:
// $(E_n-E_p)/(2\Delta{t}) = 1/\epsilon_0(-\Delta{H_z}/\Delta{x} - \sigma(E_n+E_p)/2)$
// Normalize: $E_n - E_p = -2\Delta{t}/\epsilon_0 \Delta{H_z}/\Delta{x} - \sigma/\epsilon_0 \Delta{t} (E_n + E_p)$
// Rearrange: $E_n(1 + \sigma/\epsilon_0 \Delta{t}) = E_p(1 - \sigma/\epsilon_0 \Delta{t}) - 2\Delta{t}/\epsilon_0 \Delta{H_z}/\Delta{x}$
// Then $E_n$ (i.e. the y value of $E$ after this step) is trivially solved
// ```
use consts::real::{C, C2, EPS0, HALF, ONE, TWO};
let sigma: R64 = self.mat.conductivity().into();
// XXX not obvious that bz_to_hz is sensible
let delta_hz_y = self.hz() - self.bz_to_hz(up.bz());
let ex_rhs = self.state.ex*(ONE() - sigma/EPS0()*delta_t) + TWO()*delta_t/EPS0()*delta_hz_y/feature_size;
let ex_next = ex_rhs / (ONE() + sigma/EPS0()*delta_t);
let delta_hz_x = self.hz() - self.bz_to_hz(left.bz());
let ey_rhs = self.state.ey*(ONE() - sigma/EPS0()*delta_t) - TWO()*delta_t/EPS0()*delta_hz_x/feature_size;
let ey_next = ey_rhs / (ONE() + sigma/EPS0()*delta_t);
Cell {
state: CellState {
ex: ex_next,
ey: ey_next,
hz: self.state.hz,
},
mat: self.mat,
}
}
}
#[derive(Copy, Clone, Default)]
pub struct CellState {
ex: R64,
ey: R64,
hz: R64,
}
impl CellState {
pub fn ex(&self) -> f64 {
self.ex.into()
}
pub fn ey(&self) -> f64 {
self.ey.into()
}
pub fn hz(&self) -> f64 {
self.hz.into()
}
}