From ff1342ff8a12a44f7389a74211d81a16f643bd3e Mon Sep 17 00:00:00 2001 From: colin Date: Tue, 30 Aug 2022 15:25:57 -0700 Subject: [PATCH] sim: spirv: get the `ray_propagation` test working --- crates/coremem/src/sim/spirv/mod.rs | 111 ++++++++++++---------------- 1 file changed, 48 insertions(+), 63 deletions(-) diff --git a/crates/coremem/src/sim/spirv/mod.rs b/crates/coremem/src/sim/spirv/mod.rs index abedf67..fdc480b 100644 --- a/crates/coremem/src/sim/spirv/mod.rs +++ b/crates/coremem/src/sim/spirv/mod.rs @@ -171,6 +171,14 @@ impl SpirvSim diag: SyncDiagnostics::new(), } } + + /// XXX for test + pub(crate) fn set_fields_at_index(&mut self, pos: Index, f: Fields) { + let idx = self.flat_index(pos).unwrap(); + self.e[idx] = f.e(); + self.h[idx] = f.h(); + self.m[idx] = f.m(); + } } impl SpirvSim @@ -553,77 +561,54 @@ mod test { #[test] fn ray_propagation() { - let size = Index::new(1, 1, 100); + let size = Index::new(1, 1, 1536); let feat_size = 1e-3; let mut sim = SpirvSim::, $backend>::new(size, feat_size); + let time_step = sim.timestep(); - // 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 wave = |z, t| { + let amp_e = R32::one(); + let amp_h = R32::c_inv() * R32::mu0_inv(); + // this describes a rightward-traveling wave that starts at z=512, + // with a wavelength of 512, and we only render one cycle of it. + let t_shift = t as f32 * time_step/feat_size * f32::c(); + let idx = ((z as f32 - 512.0 - t_shift) / 256.0); + if idx < 0.0 || idx > 2.0 { + return (R32::zero(), R32::zero()); + } + let phase = idx.cast::() * R32::pi(); + let a = phase.sin(); + (amp_e * a, amp_h * a) - 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)]); + }; + + // inject the wave at t=0 + for i in 0..size.z() { + let (e, h) = wave(i, 0); + sim.set_fields_at_index(Index::new(0, 0, i), Fields::new( + Vec3::new(e, R32::zero(), R32::zero()), + Vec3::new(R32::zero(), h, R32::zero()), + Vec3::default(), + )); - 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::(), f.h().cast::()) - }).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; + // advance the simulation N steps. + // at each step, the wave we injected should have been shifted by some known amount, + // and we assert that the wave is where we expect it (and with the right + // amplitude) + for t in 1..512 { + sim.step(); + let eh: Vec<_> = (0..size.z()).into_iter().map(|z| { + let f = sim.fields_at_index(Index::new(0, 0, z)); + (f.e().x().cast::(), f.h().y().cast::()) + }).collect(); + for (z, (e, h)) in eh.iter().enumerate() { + let (exp_e, exp_h) = wave(z as u32, t); + assert_float_eq!(*e, exp_e.cast::(), abs <= 1e-2, "t={}, z={} ... {:?}", t, z, eh); + assert_float_eq!(*h, exp_h.cast::(), abs <= 1e-4, "t={}, z={} ... {:?}", t, z, eh); + } } - // 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); } } }