From 432d5a032e1cb4c2bc1bf326c72ccd314196944a Mon Sep 17 00:00:00 2001 From: Colin Date: Sat, 29 May 2021 16:50:25 -0700 Subject: [PATCH] Define a `ConvFlt` type for computing the PML convolutions Not used beyond tests, yet. --- src/sim.rs | 166 +++++++++++++++++++++++++++++++++++++++++++++++++++-- 1 file changed, 161 insertions(+), 5 deletions(-) diff --git a/src/sim.rs b/src/sim.rs index 1037590..185503f 100644 --- a/src/sim.rs +++ b/src/sim.rs @@ -15,10 +15,10 @@ pub type StaticSim = SimState; #[derive(Default, Copy, Clone, PartialEq)] pub struct PmlState { - /// Auxiliary field used when stepping E - aux_e: Vec3, - /// Auxiliary field used when stepping B - aux_b: Vec3, + /// Auxiliary fields used when stepping E + aux_e: PmlAux, + /// Auxiliary fields used when stepping B + aux_b: PmlAux, } impl PmlState { @@ -27,6 +27,55 @@ impl PmlState { } } +#[derive(Default, Copy, Clone, PartialEq)] +struct PmlAux { + /// ```tex + /// $Sy^-1 \ast dFz/dy$, where F is either E or H. + /// This is used to advance the X coordinate of the opposite field + /// ``` + y_conv_dfz_dy: ConvFlt, + z_conv_dfy_dz: ConvFlt, + /// Used to advance the Y coordinate of the opposite field + z_conv_dfx_dz: ConvFlt, + x_conv_dfz_dx: ConvFlt, + /// Used to advance the Z coordinate of the opposite field + x_conv_dfy_dx: ConvFlt, + y_conv_dfx_dy: ConvFlt, +} + +#[derive(Default, Copy, Clone, PartialEq)] +struct ConvFlt(Flt); + +impl ConvFlt { + ///```tex + /// Consider $v(t) = (a \delta(t) + b exp[-c t] u(t)) \ast y(t)$ + /// This method will advance from self = $v(t)$ to self = $v(t + \Delta t)$ + /// ``` + fn step(&mut self, a: Flt, b: Flt, c: Flt, y: Flt, delta_t: Flt) -> Flt { + //```tex + // let $v(t) = a \delta(t) \ast y(t) + m(t)$ + // where $m(t) = b exp[-c t] u(t) \ast y(t)$ + // let $k(t) = b exp[-c t] u(t)$, where "k" is for "Kernel". + // $m(t) = \int_0^{t} k(T) y(t - T)dT$ + // $m(t) = \int_0^{\Delta t} k(T) y(t - T)dT + \int_{\Delta t}^t k(T) y(t - T)dT$ + // $m(t) = \int_0^{\Delta t} k(T) y(t - T)dT + \int_{\Delta t}^t b exp[-c T] y(t - T)dT$ + // $m(t) = \int_0^{\Delta t} k(T) y(t - T)dT + \int_0^{t-\Delta t} b exp[-c T] exp[-c \Delta t] y(t - T)dT$ + // $m(t) = \int_0^{\Delta t} k(T) y(t - T)dT + exp[-c \Delta t] v_prev$ + // + // Assume that $y(t)$ is constant for this time-step. + // Looking just at the integral term: + // $b y \int_0^{\Delta t} exp[-c T] dT$ + // $= b y [-1/c](exp[-c \Delta t] - 1)$ + // So + // $m(t) = exp[-c \Delta t] v_prev + -b y/c (exp[-c \Delta t] - 1)$ + //``` + let e = (-c * delta_t).exp(); + self.0 = e*self.0 - b*y/c * (e - 1.0); + self.0 + a*y + } +} + + /// Used to instruct the Sim how to advance the state at a given location. /// /// `UpdateMethod::Maxwell` means to apply Maxwell's equations as usual, with @@ -43,9 +92,9 @@ pub enum UpdateMethodMut<'a> { Pml { pml: &'a mut PmlState, pseudo_conductivity: Vec3 }, } +/// See `UpdateMethodMut` #[derive(PartialEq)] pub enum UpdateMethod<'a> { - /// `conductivity` is traditionally represented by $\sigma$ Maxwell { conductivity: Vec3 }, Pml { pml: &'a PmlState, pseudo_conductivity: Vec3 }, } @@ -876,6 +925,113 @@ mod test { use rand::rngs::StdRng; use rand::{Rng as _, SeedableRng as _}; + fn do_conv_flt_test(a: Flt, b: Flt, c: Flt, delta_t: Flt, signal: &[Flt], expected: &[Flt]) { + let mut cf = ConvFlt(0.0); + for (s, e) in signal.iter().zip(expected.iter()) { + let v = cf.step(a, b, c, *s, delta_t); + assert_float_eq!(v, e, r2nd <= 1e-6); + } + } + + #[test] + fn conv_flt_trivial() { + let signal = [2.0, 7.0, -3.0]; + // kernel: 1.0 \delta(t) + let expected = [2.0, 7.0, -3.0]; + let (a, b, c, delta_t): (Flt, Flt, Flt, Flt) = (1.0, 0.0, 1.0, 0.1); + do_conv_flt_test(a, b, c, delta_t, &signal, &expected); + } + + #[test] + fn conv_flt_trivial_scaled() { + let signal = [2.0, 7.0, -3.0]; + // kernel: 10.0 \delta(t) + let expected = [20.0, 70.0, -30.0]; + let (a, b, c, delta_t): (Flt, Flt, Flt, Flt) = (10.0, 0.0, 1.0, 0.1); + do_conv_flt_test(a, b, c, delta_t, &signal, &expected); + } + + #[test] + fn conv_flt_exp_trivial() { + let signal = [2.0, 0.0, 0.0]; + // kernel: e(-t) + // \int_0^1 e(-t) dt = [1 - exp(-1)] + let exp_neg_0 = (-0.0 as Flt).exp(); + let exp_neg_1 = (-1.0 as Flt).exp(); + let exp_neg_2 = (-2.0 as Flt).exp(); + let exp_neg_3 = (-3.0 as Flt).exp(); + let expected = [ + 2.0*(exp_neg_0 - exp_neg_1), + 2.0*(exp_neg_1 - exp_neg_2), + 2.0*(exp_neg_2 - exp_neg_3), + ]; + let (a, b, c, delta_t): (Flt, Flt, Flt, Flt) = (0.0, 1.0, 1.0, 1.0); + do_conv_flt_test(a, b, c, delta_t, &signal, &expected); + } + + #[test] + fn conv_flt_exp_time_scaled() { + let signal = [2.0, 0.0, 0.0]; + // kernel: e(-3*t) + // \int_0^0.2 e(-3*t) dt = [1 - exp(-0.6)]/3 + let exp_neg_00 = (-0.0 as Flt).exp(); + let exp_neg_06 = (-0.6 as Flt).exp(); + let exp_neg_12 = (-1.2 as Flt).exp(); + let exp_neg_18 = (-1.8 as Flt).exp(); + let expected = [ + 2.0/3.0*(exp_neg_00 - exp_neg_06), + 2.0/3.0*(exp_neg_06 - exp_neg_12), + 2.0/3.0*(exp_neg_12 - exp_neg_18), + ]; + let (a, b, c, delta_t): (Flt, Flt, Flt, Flt) = (0.0, 1.0, 3.0, 0.2); + do_conv_flt_test(a, b, c, delta_t, &signal, &expected); + } + + #[test] + fn conv_flt_exp_multi_input() { + let signal = [2.0, 7.0, -3.0]; + // kernel: e(-3*t) + // \int_0^0.2 e(-3*t) dt = [1 - exp(-0.6)]/3 + let exp_neg_00 = (-0.0 as Flt).exp(); + let exp_neg_06 = (-0.6 as Flt).exp(); + let exp_neg_12 = (-1.2 as Flt).exp(); + let exp_neg_18 = (-1.8 as Flt).exp(); + let expected = [ + 2.0/3.0*(exp_neg_00 - exp_neg_06), + 7.0/3.0*(exp_neg_00 - exp_neg_06) + 2.0/3.0*(exp_neg_06 - exp_neg_12), + -3.0/3.0*(exp_neg_00 - exp_neg_06) + + 7.0/3.0*(exp_neg_06 - exp_neg_12) + + 2.0/3.0*(exp_neg_12 - exp_neg_18), + ]; + let (a, b, c, delta_t): (Flt, Flt, Flt, Flt) = (0.0, 1.0, 3.0, 0.2); + do_conv_flt_test(a, b, c, delta_t, &signal, &expected); + } + + #[test] + fn conv_flt_all_params() { + let signal = [2.0, 7.0, -3.0]; + // kernel: 0.3 \delta(t) + 0.4 * e(-3*t) + // \int_0^0.2 e(-3*t) dt = [1 - exp(-0.6)]/3 + let exp_neg_00 = (-0.0 as Flt).exp(); + let exp_neg_06 = (-0.6 as Flt).exp(); + let exp_neg_12 = (-1.2 as Flt).exp(); + let exp_neg_18 = (-1.8 as Flt).exp(); + let expected_exp = [ + 2.0/3.0*(exp_neg_00 - exp_neg_06), + 7.0/3.0*(exp_neg_00 - exp_neg_06) + 2.0/3.0*(exp_neg_06 - exp_neg_12), + -3.0/3.0*(exp_neg_00 - exp_neg_06) + + 7.0/3.0*(exp_neg_06 - exp_neg_12) + + 2.0/3.0*(exp_neg_12 - exp_neg_18), + ]; + let expected = [ + 0.3 * signal[0] + 0.4 * expected_exp[0], + 0.3 * signal[1] + 0.4 * expected_exp[1], + 0.3 * signal[2] + 0.4 * expected_exp[2], + ]; + let (a, b, c, delta_t): (Flt, Flt, Flt, Flt) = (0.3, 0.4, 3.0, 0.2); + do_conv_flt_test(a, b, c, delta_t, &signal, &expected); + } + fn assert_vec3_eq(actual: Vec3, expected: Vec3) { assert_float_eq!(actual.x(), expected.x(), r2nd <= 1e-6); assert_float_eq!(actual.y(), expected.y(), r2nd <= 1e-6);