From 5677670d0568fd43dc9947aea820c3b23c35c82c Mon Sep 17 00:00:00 2001 From: Colin Date: Sat, 22 Aug 2020 13:50:30 -0700 Subject: [PATCH] Fix the incorrect expansion of E about T in Cell::step_e Now the coremem example can handle a conductor of any reasonable conductivity without causing numerical errors. --- src/lib.rs | 69 ++++++++++++++++++++++++++++++------------------------ 1 file changed, 38 insertions(+), 31 deletions(-) diff --git a/src/lib.rs b/src/lib.rs index c90b2ac..3a2ab92 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -31,9 +31,18 @@ pub mod consts { pub fn C() -> R64 { super::C.into() } + pub fn C2() -> R64 { + C() * C() + } pub fn MU0() -> R64 { super::MU0.into() } + pub fn ONE() -> R64 { + 1.0.into() + } + pub fn TWO() -> R64 { + 2.0.into() + } pub fn HALF() -> R64 { 0.5.into() } @@ -63,7 +72,6 @@ impl SimState { pub fn step(&mut self) { use consts::real::*; 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())); // first advance all the magnetic fields for down_y in 1..self.height() { @@ -82,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()); + 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); @@ -112,7 +120,6 @@ impl SimState { &mut self.cells[[y, x]] } - fn timestep(&self) -> R64 { self.feature_size / consts::real::C() } @@ -186,46 +193,46 @@ impl Cell { } /// 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. - /// 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) -> Self { + /// 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: // $\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$ - // 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 - // Rearrange: $dE/dt = c^2 (\nabla x B - S)$ - // Expand: $dE_x/dt = c^2 (dB_z/dy - S_x)$ (1); $dE_y/dt = c^2 (-dB_z/dx - S_y)$ (2) + // 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 - \mu_0 \sigma E)$ + // 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)$ - // Information is propagated across $1/2 \Delta x$ where $\Delta x$ = grid spacing of cells. - // Therefore $1/2 \Delta x = c \Delta t$ or $\Delta t / \Delta x = 1/(2c)$ - // Rearrange: $\Delta E_x = c^2 (\Delta B_z \Delta t / \Delta y - \Delta t S_x)$ - // Rearrange: $\Delta E_x = c (\Delta B_z/2 - c \Delta t S_x)$ - // - // 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 (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) = c^2(\Delta B_z/\Delta y - \mu_0\sigma(E_n+E_p)/2)$ + // 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: $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 + // + // 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; - // TYPO: self.ex here should actually be replaced with 0.5(self.ex + self.ex+delta_ex). - // i.e. it should be the Ex halfway through the timestep. - // 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 ex_rhs = self.ex*(ONE() - C2()*MU0()*sigma*delta_t) + TWO()*delta_t*C2()*delta_bz_y/feature_size; + let ex_next = ex_rhs / (ONE() + C2()*MU0()*sigma*delta_t); let delta_bz_x = self.bz - left.bz; - let static_ey: R64 = MU0() * self.mat.conductivity() * self.ey; - let delta_ey: R64 = C() * (-HALF() * delta_bz_x - C() * delta_t * static_ey); - + let ey_rhs = self.ey*(ONE() - C2()*MU0()*sigma*delta_t) - TWO()*delta_t*C2()*delta_bz_x/feature_size; + let ey_next = ey_rhs / (ONE() + C2()*MU0()*sigma*delta_t); + Cell { - ex: self.ex + delta_ex, - ey: self.ey + delta_ey, + ex: ex_next, + ey: ey_next, bz: self.bz, mat: self.mat, }