From c33317c80562359be09cb1199f99b22fa70d1666 Mon Sep 17 00:00:00 2001 From: Colin Date: Sat, 3 Oct 2020 15:25:30 -0700 Subject: [PATCH] Factor out a `Neighbors` to package up all 6 neighbors when updating cells This will allow more flexibility in boundary conditions. It seems to have similar perf as passing only the 3 cells within about 2%. --- src/sim.rs | 56 ++++++++++++++++++++++++++++++++++++++---------------- 1 file changed, 40 insertions(+), 16 deletions(-) diff --git a/src/sim.rs b/src/sim.rs index 64d2b9f..2bbd115 100644 --- a/src/sim.rs +++ b/src/sim.rs @@ -98,10 +98,7 @@ impl SimState { Zip::from(ndarray::indices_of(&self.cells)).and(&self.cells).par_apply_assign_into( &mut scratch_cells, |(z, y, x), cell| { - let down_cell = self.cells.get([z, y+1, x]); - let right_cell = self.cells.get([z, y, x+1]); - let in_cell = self.cells.get([z+1, y, x]); - cell.clone().step_h(right_cell, down_cell, in_cell, half_time_step.into(), feature_size.into()) + cell.clone().step_h(Neighbors::new(&self.cells, [z, y, x]), half_time_step.into(), feature_size.into()) }); trace!("step_h end"); @@ -110,10 +107,7 @@ impl SimState { Zip::from(ndarray::indices_of(&scratch_cells)).and(&scratch_cells).par_apply_assign_into( &mut self.cells, |(z, y, x), cell| { - let left_cell = scratch_cells.get([z, y, x.wrapping_sub(1)]); - let up_cell = scratch_cells.get([z, y.wrapping_sub(1), x]); - let out_cell = scratch_cells.get([z.wrapping_sub(1), y, x]); - cell.clone().step_e(left_cell, up_cell, out_cell, half_time_step.into(), feature_size.into()) + cell.clone().step_e(Neighbors::new(&scratch_cells, [z, y, x]), half_time_step.into(), feature_size.into()) }); trace!("step_e end"); @@ -247,6 +241,30 @@ pub struct Cell { mat: M, } +struct Neighbors<'a, Cell> { + left: Option<&'a Cell>, + right: Option<&'a Cell>, + up: Option<&'a Cell>, + down: Option<&'a Cell>, + out: Option<&'a Cell>, + in_: Option<&'a Cell>, +} + +impl<'a, Cell> Neighbors<'a, Cell> { + fn new(array: &'a Array3, location: [usize; 3]) -> Self { + let [z, y, x] = location; + let left = array.get([z, y, x.wrapping_sub(1)]); + let right = array.get([z, y, x + 1]); + let down = array.get([z, y.wrapping_sub(1), x]); + let up = array.get([z, y + 1, x]); + let out = array.get([z.wrapping_sub(1), y, x]); + let in_ = array.get([z + 1, y, x]); + Self { + left, right, up, down, out, in_ + } + } +} + impl Cell { pub fn ex(&self) -> Flt { self.state.e().x() @@ -297,7 +315,7 @@ impl Cell { self.state.h = self.mat.step_h(&self.state, Vec3::new(0.0, 0.0, delta_bz)); } - fn step_h(mut self, right: Option<&Self>, down: Option<&Self>, in_: Option<&Self>, delta_t: Real, feature_size: Real) -> Self { + fn step_h(mut self, neighbors: Neighbors<'_, Self>, delta_t: Real, feature_size: Real) -> Self { // ```tex // Maxwell-Faraday equation: $\nabla x E = -dB/dt$ // @@ -313,21 +331,21 @@ impl Cell { use consts::real::ZERO; let inv_feature_size = Real::from_inner(1.0) / feature_size; - let (delta_ey_x, delta_ez_x) = match right { + let (delta_ey_x, delta_ez_x) = match neighbors.right { None => (ZERO(), ZERO()), Some(right) => ( inv_feature_size * (right.ey() - self.ey()), inv_feature_size * (right.ez() - self.ez()), ) }; - let (delta_ex_y, delta_ez_y) = match down { + let (delta_ex_y, delta_ez_y) = match neighbors.down { None => (ZERO(), ZERO()), Some(down) => ( inv_feature_size * (down.ex() - self.ex()), inv_feature_size * (down.ez() - self.ez()), ) }; - let (delta_ex_z, delta_ey_z) = match in_ { + let (delta_ex_z, delta_ey_z) = match neighbors.in_ { None => (ZERO(), ZERO()), Some(in_) => ( inv_feature_size * (in_.ex() - self.ex()), @@ -349,7 +367,7 @@ impl Cell { /// delta_t = timestep covered by this function. i.e. it should be half the timestep of the simulation /// feature_size = how many units apart is the center of each adjacent cell on the grid. - fn step_e(self, left: Option<&Self>, up: Option<&Self>, out: Option<&Self>, delta_t: Real, feature_size: Real) -> Self { + fn step_e(self, neighbors: Neighbors<'_, Self>, delta_t: Real, feature_size: Real) -> Self { // ```tex // Ampere's circuital law with Maxwell's addition, in SI units ("macroscopic version"): // $\nabla x H = J_f + dD/dt$ where $J_f$ = current density = $\sigma E$, $\sigma$ being a material parameter ("conductivity") @@ -371,6 +389,12 @@ impl Cell { // | d/dx d/dy d/dz | = | dHx/dz - dHz/dx | // | | | | // | Hx Hy Hz | [ dHy/dx - dHx/dy ] + // + // Boundary conditions (TODO): + // For something on the yz plane, E(x=0, t+1) should be E(x=1, t). i.e. the wave travels + // delta_x = 1 per timestep. Roughly, we want to take the line from E(x=0, t), E(x=1, t) + // and continue that to E(x=-1, t) ? + // // ``` use consts::real::{EPS0, TWO, ZERO}; @@ -378,21 +402,21 @@ impl Cell { let sigma: Vec3 = self.mat.conductivity(); let inv_feature_size = Real::from_inner(1.0) / feature_size; - let (delta_hy_x, delta_hz_x) = match left { + let (delta_hy_x, delta_hz_x) = match neighbors.left { None => (ZERO(), ZERO()), Some(left) => ( inv_feature_size * (self.hy() - left.hy()), inv_feature_size * (self.hz() - left.hz()), ) }; - let (delta_hx_y, delta_hz_y) = match up { + let (delta_hx_y, delta_hz_y) = match neighbors.up { None => (ZERO(), ZERO()), Some(up) => ( inv_feature_size * (self.hx() - up.hx()), inv_feature_size * (self.hz() - up.hz()), ) }; - let (delta_hx_z, delta_hy_z) = match out { + let (delta_hx_z, delta_hy_z) = match neighbors.out { None => (ZERO(), ZERO()), Some(out) => ( inv_feature_size * (self.hx() - out.hx()),