cross: step: annotate step_e with a maths derivation

this derivation comes from the legacy cpu code.
it looks like the spirv implementation actually
pursued a conceptually simpler approach.
i'll try to reconcile these shortly.
This commit is contained in:
2022-08-25 22:20:03 -07:00
parent 7f40b3ccd5
commit bbd31cc7be

View File

@@ -65,19 +65,41 @@ impl<'a, R: Real, M: Material<R>> StepEContext<'a, R, M> {
/// return the e field for this cell at `t=t_{now} + \delta t`
pub fn step_e(self) -> Vec3<R> {
// ```tex
// Ampere's circuital law with Maxwell's addition, in SI units ("macroscopic version"):
// $\nabla x H = J_f + dD/dt$ (1)
// 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$.
//
// substitute $D = \epsilon_0 E$ into (1):
// $\nabla x H = J_f + \epsilon_0 dE/dt$
// expand with $J_f = \sigma E$:
// $\nabla x H = \sigma E + \epsilon_0 dE/dt$
// rearrange:
// $dE/dt = 1/\epsilon_0 (\nabla x H - \sigma E)$ (2)
//
// let $E_p$ be $E$ at $T - \Delta\!t$ and $E_n$ be $E$ at $T+\Delta\!t$.
// apply these substitutions into a linear expansion of (2):
// $(E_n - E_p)/(2\Delta\!t) = 1/\epsilon_0 (\nabla x H - \sigma (E_n + E_p)/2)$
// normalize:
// $E_n - E_p = 2\Delta\!t/\epsilon_0 (\nabla x H - \sigma (E_n + E_p)/2)$
// expand:
// $E_n - E_p = 2\Delta\!t/\epsilon_0 \nabla x H - \sigma \Delta\!t/\epsilon_0 (E_n + E_p)$
// rearrange:
// $E_n(1 + \sigma \Delta\!t/\epsilon_0) = E_p(1 - \sigma \Delta\!t/\epsilon_0) + 2\Delta\!t/\epsilon_0 \nabla x H$
// Then $E_n$ (i.e. the E vector after this step) is trivially solved
// ```
let twice_eps0 = R::twice_eps0();
let deltas = self.in_h.delta_h();
// \nabla x H
let nabla_h = deltas.nabla() * self.inv_feature_size;
// $\nabla x H = \epsilon_0 dE/dt + \sigma E$
// no-conductivity version:
// let delta_e = nabla_h * (self.time_step * EPS0_INV);
let sigma = self.mat.conductivity();
let e_prev = self.in_e;
let delta_e = (nabla_h - e_prev.elem_mul(sigma)).elem_div(
sigma*self.time_step + Vec3::uniform(twice_eps0)
)*(R::two()*self.time_step);
// println!("spirv-step_e delta_e: {:?}", delta_e);
e_prev + delta_e + self.stim_e
}
}