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

Fix the incorrect expansion of E about T in Cell::step_e

Browse Source 
Now the coremem example can handle a conductor of any reasonable
conductivity without causing numerical errors.
This commit is contained in:
Colin
2020-08-22 13:50:30 -07:00
parent dff8166906
commit 5677670d05
1 changed files with 38 additions and 31 deletions
Download Patch File Download Diff File Expand all files Collapse all files

69
src/lib.rs
View File

@@ -31,9 +31,18 @@ pub mod consts {
pub fn C() -> R64 { pub fn C() -> R64 {
super::C.into() super::C.into()
} }
pub fn C2() -> R64 {
C() * C()
}
pub fn MU0() -> R64 { pub fn MU0() -> R64 {
super::MU0.into() super::MU0.into()
} }
pub fn ONE() -> R64 {
1.0.into()
}
pub fn TWO() -> R64 {
2.0.into()
}
pub fn HALF() -> R64 { pub fn HALF() -> R64 {
0.5.into() 0.5.into()
} }
@@ -63,7 +72,6 @@ impl SimState {
pub fn step(&mut self) { pub fn step(&mut self) {
use consts::real::*; use consts::real::*;
let half_time_step = HALF() * self.timestep(); let half_time_step = HALF() * self.timestep();
let half_time_step = 0.0005 / consts::C;
let mut working_cells = Array2::default((self.height(), self.width())); let mut working_cells = Array2::default((self.height(), self.width()));
// first advance all the magnetic fields // first advance all the magnetic fields
for down_y in 1..self.height() { for down_y in 1..self.height() {
@@ -82,7 +90,7 @@ impl SimState {
let cell = self.get(x, y); let cell = self.get(x, y);
let left_cell = self.get(x-1, y); let left_cell = self.get(x-1, y);
let up_cell = self.get(x, y-1); let up_cell = self.get(x, y-1);
working_cells[[y, x]] = cell.step_e(left_cell, up_cell, half_time_step.into()); working_cells[[y, x]] = cell.step_e(left_cell, up_cell, half_time_step.into(), self.feature_size.into());
} }
} }
std::mem::swap(&mut working_cells, &mut self.cells); std::mem::swap(&mut working_cells, &mut self.cells);
@@ -112,7 +120,6 @@ impl SimState {
&mut self.cells[[y, x]] &mut self.cells[[y, x]]
} }
fn timestep(&self) -> R64 { fn timestep(&self) -> R64 {
self.feature_size / consts::real::C() self.feature_size / consts::real::C()
} }
@@ -186,46 +193,46 @@ impl<M: Material + Clone> Cell<M> {
} }
/// delta_t = timestep covered by this function. i.e. it should be half the timestep of the simulation /// delta_t = timestep covered by this function. i.e. it should be half the timestep of the simulation
/// since the simulation spends half a timestep in step_b and the other half in step_e. /// feature_size = how many units apart is the center of each adjacent cell on the grid.
/// delta_x and delta_y are derived from delta_t (so, make sure delta_t is constant across all calls if the grid spacing is also constant!) 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) -> Self {
// ```tex // ```tex
// Ampere's circuital law with Maxwell's addition, in SI units: // 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 // $\nabla x B = \mu_0 (J + \epsilon_0 dE/dt)$ where J = current density = $\sigma E$, $\sigma$ being a material parameter
// Expand: $\nabla x B = \mu_0 \sigma E + \mu_0 \epsilon_0 dE/dt$ // Expand: $\nabla x B = \mu_0 \sigma E + \mu_0 \epsilon_0 dE/dt$
// Substitute: $\nabla x B = S + 1/c^2 dE/dt$ where $c = 1/\sqrt{\mu_0 \epsilon_0}$ is the speed of light, and $S = \mu_0 \sigma E$ for convenience // Substitute: $\nabla x B = \mu_0 \sigma E + 1/c^2 dE/dt$ where $c = 1/\sqrt{\mu_0 \epsilon_0}$ is the speed of light
// Rearrange: $dE/dt = c^2 (\nabla x B - S)$ // Rearrange: $dE/dt = c^2 (\nabla x B - \mu_0 \sigma E)$
// Expand: $dE_x/dt = c^2 (dB_z/dy - S_x)$ (1); $dE_y/dt = c^2 (-dB_z/dx - S_y)$ (2) // Expand: $dE_x/dt = c^2 (dB_z/dy - \mu_0 \sigma E_x)$ (1); $dE_y/dt = c^2 (-dB_z/dx - \mu_0 \sigma E_y)$ (2)
// //
// Discretize (1): $(\Delta E_x)/(\Delta t) = c^2 (\Delta B_z / \Delta y - S_x)$ // Consider (1): let $E_p$ be $E_x$ at $T-\Delta t$ and $E_n$ be $E_x$ at $T+\Delta t$.
// Information is propagated across $1/2 \Delta x$ where $\Delta x$ = grid spacing of cells. // Linear expansion about $t=T$, and discretized:
// Therefore $1/2 \Delta x = c \Delta t$ or $\Delta t / \Delta x = 1/(2c)$ // $(E_n-E_p)/(2\Delta t) = c^2(\Delta B_z/\Delta y - \mu_0\sigma(E_n+E_p)/2)$
// Rearrange: $\Delta E_x = c^2 (\Delta B_z \Delta t / \Delta y - \Delta t S_x)$ // Normalize: $E_n - E_p = 2\Delta{t} c^2 \Delta{B_z}/\Delta{y} - c^2 \mu_0 \sigma \Delta{t} (E_n + E_p)$
// Rearrange: $\Delta E_x = c (\Delta B_z/2 - c \Delta t S_x)$ // Rearrange: $E_n(1 + c^2 \mu_0 \sigma \Delta{t}) = E_p(1 - c^2 \mu_0 \sigma \Delta{t}) + 2\Delta{t} c^2 \Delta{B_z}/\Delta{y}$
// // Then $E_n$ (i.e. the x value of $E$ after this step) is trivially solved
// Discretize (2): $(\Delta E_y)/(\Delta t) = c^2 (-\Delta B_z / \Delta x - S_y)$ //
// Rearrange: $\Delta E_y = c (-\Delta B_z / 2 - c \Delta_t S_y)$ // 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) = c^2(-\Delta B_z/\Delta x - \mu_0\sigma(E_n+E_p)/2)$
// Normalize: $E_n - E_p = -2\Delta{t} c^2 \Delta{B_z}/\Delta{x} - c^2 \mu_0 \sigma \Delta{t} (E_n + E_p)$
// Rearrange: $E_n(1 + c^2 \mu_0 \sigma \Delta{t}) = E_p(1 - c^2 \mu_0 \sigma \Delta{t}) - 2\Delta{t} c^2 \Delta{B_z}/\Delta{x}$
// Then $E_n$ (i.e. the y value of $E$ after this step) is trivially solved
// ``` // ```
use consts::real::{C, HALF, MU0}; use consts::real::{C, C2, HALF, MU0, ONE, TWO};
let sigma = self.mat.conductivity();
let delta_bz_y = self.bz - up.bz; let delta_bz_y = self.bz - up.bz;
// TYPO: self.ex here should actually be replaced with 0.5(self.ex + self.ex+delta_ex). let ex_rhs = self.ex*(ONE() - C2()*MU0()*sigma*delta_t) + TWO()*delta_t*C2()*delta_bz_y/feature_size;
// i.e. it should be the Ex halfway through the timestep. let ex_next = ex_rhs / (ONE() + C2()*MU0()*sigma*delta_t);
// BUT: meep defines metals as a medium with epsilon=-\inf. No mention of
// conductivity/sigma. Maybe that route is simpler (if equivalent?)
// Note that metals don't really seem to have 'bound' current: just 'free' current:
// https://physics.stackexchange.com/questions/227014/are-the-conducting-electrons-in-a-metal-counted-as-free-or-bound-charges
let static_ex: R64 = MU0() * self.mat.conductivity() * self.ex;
let delta_ex: R64 = C() * (HALF() * delta_bz_y - C() * delta_t * static_ex);
let delta_bz_x = self.bz - left.bz; let delta_bz_x = self.bz - left.bz;
let static_ey: R64 = MU0() * self.mat.conductivity() * self.ey; let ey_rhs = self.ey*(ONE() - C2()*MU0()*sigma*delta_t) - TWO()*delta_t*C2()*delta_bz_x/feature_size;
let delta_ey: R64 = C() * (-HALF() * delta_bz_x - C() * delta_t * static_ey); let ey_next = ey_rhs / (ONE() + C2()*MU0()*sigma*delta_t);
Cell { Cell {
ex: self.ex + delta_ex, ex: ex_next,
ey: self.ey + delta_ey, ey: ey_next,
bz: self.bz, bz: self.bz,
mat: self.mat, mat: self.mat,
} }