add a chaos UPML test (broken)

Cargo.toml

@@ -27,6 +27,7 @@ ndarray = { version = "0.13", features = ["rayon", "serde"] }
 num = "0.3"
 piecewise-linear = "0.1"
 plotly = { version = "0.6", features = ["kaleido", "plotly_ndarray"], path = "../plotly/plotly" }
+rand = "0.8"
 rayon = "1.5"
 serde = "1.0"
 structopt = "0.3"

src/sim.rs

@@ -798,6 +798,8 @@ mod test {
     use super::*;
     use crate::geom::WorldRegion;
     use crate::mat::Conductor;
+    use rand::rngs::StdRng;
+    use rand::{Rng as _, SeedableRng as _};
     use std::sync::Mutex;
 
     fn energy(s: &dyn GenericSim) -> Flt {
@@ -1009,8 +1011,6 @@ mod test {
         upml_test(&mut state, Vec3::unit_z());
     }
 
-    // TODO:
-    //   3d UPML chaos testing (set the E/B field to completely random state)
     #[test]
     fn upml_boundary_3d_multidirectional() {
         let end = 40;
@@ -1041,4 +1041,64 @@ mod test {
         upml_test(&mut state, Vec3::unit_y());
         upml_test(&mut state, Vec3::unit_z());
     }
+
+    #[test]
+    fn upml_boundary_3d_chaotic_field() {
+        let end = 40;
+        let mut state = StaticSim::new(Index((end+1, end+1, end+1).into()), 1e-6);
+        for inset in 0..end/4 {
+            let f = 1e6 * (end/4 - inset) as Flt;
+            let stretch_x = Vec3::new(f*f, 0.0, 0.0);
+            let stretch_y = Vec3::new(0.0, f*f, 0.0);
+            let stretch_z = Vec3::new(0.0, 0.0, f*f);
+            for in_plane_a in 0..=end {
+                for in_plane_b in 0..=end {
+                    state.get_mut(Index::new(inset, in_plane_a, in_plane_b))
+                        .mat_mut().coord_stretch += stretch_x;
+                    state.get_mut(Index::new(end-inset, in_plane_a, in_plane_b))
+                        .mat_mut().coord_stretch += stretch_x;
+                    state.get_mut(Index::new(in_plane_a, inset, in_plane_b))
+                        .mat_mut().coord_stretch += stretch_y;
+                    state.get_mut(Index::new(in_plane_a, end-inset, in_plane_b))
+                        .mat_mut().coord_stretch += stretch_y;
+                    state.get_mut(Index::new(in_plane_a, in_plane_b, inset))
+                        .mat_mut().coord_stretch += stretch_z;
+                    state.get_mut(Index::new(in_plane_a, in_plane_b, end-inset))
+                        .mat_mut().coord_stretch += stretch_z;
+                }
+            }
+        }
+        // psuedo-randomly perturb the field:
+        for t in 0..100 {
+            // re-seeding the rng allows each coordinate to make the same choices
+            // across all time steps
+            let mut rng = StdRng::seed_from_u64(0x19b562c7_fa03be59);
+            let angle = ((t as Flt)/100.0*2.0*consts::PI);
+            let gate = 0.5*(1.0 - angle.cos());
+            for z in 0..=end {
+                for y in 0..=end {
+                    for x in 0..=end {
+                        let omega = rng.gen_range(0.01..0.25) * 2.0*consts::PI;
+                        let sig = (omega * t as Flt).sin();
+                        let dir = Vec3::new(
+                            rng.gen_range(-1.0..1.0),
+                            rng.gen_range(-1.0..1.0),
+                            rng.gen_range(-1.0..1.0),
+                        );
+                        let excitation = dir * sig * gate;
+                        let dyn_state = &mut state as &mut dyn GenericSim;
+                        dyn_state.impulse_e(Index::new(x, y, z), excitation);
+                    }
+                }
+            }
+        }
+        // let the energy drain
+        let energy_0 = energy(&state);
+        for _ in 0..2000 {
+            state.step();
+        }
+        let energy_1 = energy(&state);
+        // PML should absorb all energy
+        assert_float_eq!(energy_1/energy_0, 0.0, abs <= 1e-5);
+    }
 }