spirv: test: prototype a ray_propagation test

it needs some more work to be useful

crates/coremem/src/sim/spirv/mod.rs

@@ -344,12 +344,17 @@ mod test {
             mod $shortname {
                 use super::*;
 
-                use crate::geom::WorldRegion;
-                use crate::stim::{NoopStimulus, RngStimulus};
+                use crate::geom::{Coord as _, WorldRegion};
+                use crate::sim::Cube;
+                use crate::stim::{DynStimuli, FieldMags, NoopStimulus, RegionGated, RngStimulus};
                 use coremem_cross::mat::{AnisomorphicConductor, IsomorphicConductor};
                 use float_eq::assert_float_eq;
                 use more_asserts::{assert_gt, assert_lt};
 
+                ////////// SMOKE/SANITY TESTS //////////
+                // these test that everything is roughly wired together right, but don't ensure
+                // that the simulations will particularly reflect reality.
+
                 #[test]
                 fn accessors() {
                     let size = Index::new(15, 17, 19);
@@ -542,6 +547,84 @@ mod test {
                 fn mh_ferromagnet_smoke_100_steps_larger() {
                     mh_ferromagnet_smoke_test(0x1234, 100, Index::new(328, 252, 160));
                 }
+
+                ////////// APPLICATION TESTS //////////
+                // these test that the simulation accurately models *specific* applications
+
+                #[test]
+                fn ray_propagation() {
+                    let size = Index::new(1, 1, 100);
+                    let feat_size = 1e-3;
+                    let mut sim = SpirvSim::<R32, FullyGenericMaterial<R32>, $backend>::new(size, feat_size);
+
+                    // we'll inject this stimulus at a single point
+                    let region = Cube::new_centered(
+                        Index::new(0, 0, 50).to_meters(feat_size),
+                        Index::new(1, 1, 1).to_meters(feat_size),
+                    );
+
+                    let duration_frames = 40;
+                    let duration = 40.0 * sim.timestep();
+                    let inv_ts = sim.timestep().inv().to_r32();
+                    // let f = inv_ts * 0.05.to_r32();
+                    // see https://en.wikipedia.org/wiki/Electromagnetic_radiation
+                    let amp_e = Sinusoid::from_wavelength(duration.to_r32())
+                        .scaled(FieldMags::new_e(inv_ts));
+                    // since H is shifted a half grid-space, we delay it by one half timestep
+                    let amp_h = Sinusoid::from_wavelength(duration.to_r32())
+                        // .shifted(sim.timestep().to_r32() * -R32::half())
+                        .scaled(FieldMags::new_h(inv_ts * R32::c_inv() * R32::mu0_inv()));
+                    // h will be a cosine wave
+                    // let amp_h = Sinusoid::new(1.0, f).shifted(0.25 * duration);
+                    let stim_e = ModulatedVectorField::new(
+                        RegionGated::new(region.clone(), stim::Fields::new_e(Vec3::new(R32::one(), R32::zero(), R32::zero()))),
+                        amp_e,
+                    );
+                    let stim_h = ModulatedVectorField::new(
+                        RegionGated::new(region.clone(), stim::Fields::new_e(Vec3::new(R32::zero(), R32::one(), R32::zero()))),
+                        amp_h,
+                    );
+                    // let stim_h = ModulatedVectorField::new(
+                    //     amp_h,
+                    //     RegionGated::new(region.clone(), stim::Fields::new_e(Vec3::new(0.0, 0.0, 1.0))),
+                    // );
+                    let stim = DynStimuli::from_vec(vec![Box::new(stim_e), Box::new(stim_h)]);
+
+                    for _ in 0..duration_frames {
+                        sim.step_multiple(1, &stim);
+                    }
+
+                    // we should now have one wavelength of the wave propagated along the z axis.
+                    // TODO: assert correct sim state
+                    // this looks about right, just need to scale it and assert against expected
+                    // values
+                    let mut eh: Vec<_> = (0..100).into_iter().map(|z| {
+                        let f = sim.fields_at_index(Index::new(0, 0, z));
+                        (f.e().cast::<f32>(), f.h().cast::<f32>())
+                    }).collect();
+
+                    // normalize so peaks are 1.0
+                    let amp_e = eh.iter().map(|(e, _h)| e.x().abs())
+                        .fold(0f32, |a, b| a.max(b));
+                    let amp_h = eh.iter().map(|(_e, h)| h.y().abs())
+                        .fold(0f32, |a, b| a.max(b));
+                    for i in &mut eh {
+                        i.0 /= amp_e;
+                        i.1 /= amp_h;
+                    }
+                    // unpack the Ex and Hy components
+                    let mut eh_scalars = Vec::new();
+                    for i in eh {
+                        // Ey, Ez == 0
+                        assert_float_eq!(i.0.y(), 0.0, abs <= 1e-6);
+                        assert_float_eq!(i.0.z(), 0.0, abs <= 1e-6);
+                        // Hx, Hz == 0
+                        assert_float_eq!(i.1.x(), 0.0, abs <= 1e-6);
+                        assert_float_eq!(i.1.z(), 0.0, abs <= 1e-6);
+                        eh_scalars.push((i.0.x(), i.1.y()));
+                    }
+                    panic!("{:?}", eh_scalars);
+                }
             }
         }
     }