Define a ConvFlt type for computing the PML convolutions
Not used beyond tests, yet.
This commit is contained in:
166
src/sim.rs
166
src/sim.rs
@@ -15,10 +15,10 @@ pub type StaticSim = SimState<mat::Static>;
|
||||
|
||||
#[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);
|
||||
|
||||
Reference in New Issue
Block a user