colin/fdtd-coremem
1
0
Fork 0
Code Issues Pull Requests Projects Releases Wiki Activity

sim: spirv: get the ray_propagation test working

Browse Source
This commit is contained in:
colin
2022-08-30 15:25:57 -07:00
parent c050b0406f
commit ff1342ff8a
1 changed files with 48 additions and 63 deletions
Download Patch File Download Diff File Expand all files Collapse all files

111
crates/coremem/src/sim/spirv/mod.rs
View File

@@ -171,6 +171,14 @@ impl<R: Real, M: Default, B> SpirvSim<R, M, B>
diag: SyncDiagnostics::new(),
}
}
/// XXX for test
pub(crate) fn set_fields_at_index(&mut self, pos: Index, f: Fields<R>) {
let idx = self.flat_index(pos).unwrap();
self.e[idx] = f.e();
self.h[idx] = f.h();
self.m[idx] = f.m();
}
}
impl<R: Real, M: Default, B: Default> SpirvSim<R, M, B>
@@ -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::<R32, FullyGenericMaterial<R32>, $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>() * 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::<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;
// 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::<f32>(), f.h().y().cast::<f32>())
}).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::<f32>(), abs <= 1e-2, "t={}, z={} ... {:?}", t, z, eh);
assert_float_eq!(*h, exp_h.cast::<f32>(), 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);
}
}
}