From cb0c899194e69de77faa91cc5bc306bf5b82e438 Mon Sep 17 00:00:00 2001 From: Colin Date: Sat, 29 May 2021 13:00:29 -0700 Subject: [PATCH] Remove PML There are still a few references to it, but it doesn't really exist in the simulation proper. The existing implementation is confusing and ineffective (broken?). Upcoming patches will work to replace it with a more workable implementation. Also moves the Material::is_vacuum method to be on a helper trait. That's just a point of hygeine. Meant to split it out to a separate patch but messed up (-: --- src/bin/pml_tuning.rs | 2 +- src/mat.rs | 24 +++--- src/render.rs | 2 +- src/sim.rs | 166 ++++++++---------------------------------- 4 files changed, 42 insertions(+), 152 deletions(-) diff --git a/src/bin/pml_tuning.rs b/src/bin/pml_tuning.rs index d183397..0b7dccd 100644 --- a/src/bin/pml_tuning.rs +++ b/src/bin/pml_tuning.rs @@ -108,7 +108,7 @@ fn state_for_pml(size: Index, boundary: Index, feat_size: Flt, sc_coeff: Flt, co let coord_stretch = b * (sc_coeff / timestep); let conductivity = Vec3::unit() * (b.mag() * cond_coeff / timestep); Static { - coord_stretch, + // TODO PML coord_stretch, conductivity, ..Default::default() } diff --git a/src/mat.rs b/src/mat.rs index d9e9cf3..80d2c53 100644 --- a/src/mat.rs +++ b/src/mat.rs @@ -10,9 +10,6 @@ use std::cmp::Ordering; #[enum_dispatch] pub trait Material { - fn is_vacuum(&self) -> bool { - self.conductivity() == Default::default() - } /// Return \sigma, the electrical conductivity. /// For a vacuum, this is zero. For a perfect conductor, \inf. fn conductivity(&self) -> Vec3 { @@ -22,12 +19,6 @@ pub trait Material { fn m(&self) -> Vec3 { Vec3::zero() } - /// Return \gamma, which corresponds to the slope of the imaginary part of the x/y/z axis - /// wherever this material is placed. This is intended to provide a way to implement Perfectly - /// Matched Layers: "materials" which dissipate energy without reflecting it. - fn coord_stretch(&self) -> Vec3 { - Vec3::zero() - } /// Return the new h, and optionally change any internal state (e.g. magnetization). fn step_h(&mut self, context: &CellState, delta_b: Vec3) -> Vec3 { // B = mu0*(H+M), and in a default material M=0. @@ -35,13 +26,22 @@ pub trait Material { } } +pub trait MaterialExt { + fn is_vacuum(&self) -> bool; +} + +impl MaterialExt for M { + fn is_vacuum(&self) -> bool { + self.conductivity() == Default::default() + } +} + /// Capable of capturing all field-related information about a material at any /// snapshot moment-in-time. Useful for serializing state. #[derive(Clone, Default, Serialize, Deserialize)] pub struct Static { pub conductivity: Vec3, pub m: Vec3, - pub coord_stretch: Vec3, } impl Static { @@ -49,7 +49,6 @@ impl Static { Self { conductivity: m.conductivity(), m: m.m(), - coord_stretch: m.coord_stretch(), } } } @@ -61,9 +60,6 @@ impl Material for Static { fn m(&self) -> Vec3 { self.m } - fn coord_stretch(&self) -> Vec3 { - self.coord_stretch - } } impl From for Static diff --git a/src/render.rs b/src/render.rs index 21561a1..5dc7372 100644 --- a/src/render.rs +++ b/src/render.rs @@ -1,5 +1,5 @@ use crate::geom::{Index, Meters, Vec2, Vec3, Vec3u}; -use crate::{flt::{Flt, Real}, Material as _}; +use crate::{flt::{Flt, Real}, Material as _, MaterialExt as _}; use crate::mat; use crate::sim::{Cell, GenericSim, StaticSim}; use crate::meas::AbstractMeasurement; diff --git a/src/sim.rs b/src/sim.rs index 5d52b8a..8947e60 100644 --- a/src/sim.rs +++ b/src/sim.rs @@ -516,59 +516,6 @@ impl Cell { self.h().z() } - fn int_e(&self) -> Vec3 { - self.state.e_integral - } - fn int_ex(&self) -> Flt { - self.int_e().x() - } - fn int_ey(&self) -> Flt { - self.int_e().y() - } - fn int_ez(&self) -> Flt { - self.int_e().z() - } - - fn int_b(&self) -> Vec3 { - self.state.b_integral - } - fn int_bx(&self) -> Flt { - self.int_b().x() - } - fn int_by(&self) -> Flt { - self.int_b().y() - } - fn int_bz(&self) -> Flt { - self.int_b().z() - } - - fn int_h(&self) -> Vec3 { - self.state.h_integral - } - fn int_hx(&self) -> Flt { - self.int_h().x() - } - fn int_hy(&self) -> Flt { - self.int_h().y() - } - fn int_hz(&self) -> Flt { - self.int_h().z() - } - - - fn int_int_e(&self) -> Vec3 { - self.state.e_integral_integral - } - fn int_int_ex(&self) -> Flt { - self.int_int_e().x() - } - fn int_int_ey(&self) -> Flt { - self.int_int_e().y() - } - fn int_int_ez(&self) -> Flt { - self.int_int_e().z() - } - pub fn mat(&self) -> &M { &self.mat } @@ -717,59 +664,43 @@ impl Cell { let twice_delta_t = TWO() * delta_t; let half_delta_t = HALF() * delta_t; - // gamma is how a material specifies its stretched coordinates; psi, g and h are used - // internally to implement Maxwell's equations for such stretched coordinates. - let gamma = self.mat.coord_stretch(); - - let gamma_sum = Vec3::new( - gamma.y() + gamma.z(), - gamma.x() + gamma.z(), - gamma.x() + gamma.y(), - ); - let gamma_prod = Vec3::new( - gamma.y() * gamma.z(), - gamma.x() * gamma.z(), - gamma.x() * gamma.y(), - ); - let h_prev = self.h(); - let h_int_prev = self.int_h(); let b_prev = self.b(); - let b_int_prev = self.int_b(); let (delta_ey_x, delta_ez_x) = match (neighbors.right, neighbors.left) { (Some(right), _) => ( - inv_feature_size * ((right.ey() - self.ey()) + gamma.y() * (right.int_ey() - self.int_ey())), - inv_feature_size * ((right.ez() - self.ez()) + gamma.z() * (right.int_ez() - self.int_ez())), + inv_feature_size * (right.ey() - self.ey()), + inv_feature_size * (right.ez() - self.ez()), ), _ => (ZERO(), ZERO()), }; let (delta_ex_y, delta_ez_y) = match (neighbors.down, neighbors.up) { (Some(down), _) => ( - inv_feature_size * ((down.ex() - self.ex()) + gamma.x() * (down.int_ex() - self.int_ex())), - inv_feature_size * ((down.ez() - self.ez()) + gamma.z() * (down.int_ez() - self.int_ez())), + inv_feature_size * (down.ex() - self.ex()), + inv_feature_size * (down.ez() - self.ez()), ), _ => (ZERO(), ZERO()), }; let (delta_ex_z, delta_ey_z) = match (neighbors.in_, neighbors.out) { (Some(in_), _) => ( - inv_feature_size * ((in_.ex() - self.ex()) + gamma.x() * (in_.int_ex() - self.int_ex())), - inv_feature_size * ((in_.ey() - self.ey()) + gamma.y() * (in_.int_ey() - self.int_ey())), + inv_feature_size * (in_.ex() - self.ex()), + inv_feature_size * (in_.ey() - self.ey()), ), _ => (ZERO(), ZERO()), }; - // \nabla x (E + diag(\gamma) \int E) + // \nabla x E let nabla_e = Vec3::from((delta_ez_y - delta_ey_z, delta_ex_z - delta_ez_x, delta_ey_x - delta_ex_y)); + // From: $\nabla x E = -dB/dt$ let db_dt = solve_step_diff_eq( - -nabla_e, - Vec3::unit(), // dB/dt coeff - gamma_sum, // B coeff - gamma_prod, // \int B coeff + nabla_e, + -Vec3::unit(), // dB/dt coeff + Vec3::zero(), // B coeff + Vec3::zero(), // \int B coeff Vec3::zero(), // \int \int B coeff delta_t, b_prev, - b_int_prev, + Vec3::zero(), // b_int_prev, Vec3::zero(), // b_int_int_prev ); @@ -779,8 +710,6 @@ impl Cell { Cell { state: CellState { h: h.into(), - b_integral: b_int_prev + (b_prev + db_dt*delta_t)*twice_delta_t, - h_integral: h_int_prev + (h_prev + h) * delta_t, ..self.state }, mat: self.mat, @@ -824,79 +753,50 @@ impl Cell { let sigma = self.mat.conductivity(); - let gamma = self.mat.coord_stretch(); - - let gamma_sum = Vec3::new( - gamma.y() + gamma.z(), - gamma.x() + gamma.z(), - gamma.x() + gamma.y(), - ); - let gamma_prod = Vec3::new( - gamma.y() * gamma.z(), - gamma.x() * gamma.z(), - gamma.x() * gamma.y(), - ); let e_prev = self.e(); - let e_int_prev = self.int_e(); - let e_int_int_prev = self.int_int_e(); let (delta_hy_x, delta_hz_x) = match (neighbors.left, neighbors.right) { (Some(left), _) => ( - inv_feature_size * ((self.hy() - left.hy()) + gamma.y() * (self.int_hy() - left.int_hy())), - inv_feature_size * ((self.hz() - left.hz()) + gamma.z() * (self.int_hz() - left.int_hz())), + inv_feature_size * (self.hy() - left.hy()), + inv_feature_size * (self.hz() - left.hz()), ), _ => (ZERO(), ZERO()), }; let (delta_hx_y, delta_hz_y) = match (neighbors.up, neighbors.down) { (Some(up), _) => ( - inv_feature_size * ((self.hx() - up.hx()) + gamma.x() * (self.int_hx() - up.int_hx())), - inv_feature_size * ((self.hz() - up.hz()) + gamma.z() * (self.int_hz() - up.int_hz())), + inv_feature_size * (self.hx() - up.hx()), + inv_feature_size * (self.hz() - up.hz()), ), _ => (ZERO(), ZERO()), }; let (delta_hx_z, delta_hy_z) = match (neighbors.out, neighbors.in_) { (Some(out), _) => ( - inv_feature_size * ((self.hx() - out.hx()) + gamma.x() * (self.int_hx() - out.int_hx())), - inv_feature_size * ((self.hy() - out.hy()) + gamma.y() * (self.int_hy() - out.int_hy())), + inv_feature_size * (self.hx() - out.hx()), + inv_feature_size * (self.hy() - out.hy()), ), _ => (ZERO(), ZERO()), }; - // \nabla x (H + diag(\gamma) \int H) + // \nabla x H let nabla_h = Vec3::from((delta_hz_y - delta_hy_z, delta_hx_z - delta_hz_x, delta_hy_x - delta_hx_y)); + + // From: $\nabla x H = \epsilon_0 dE/dt + \sigma E$ let de_dt = solve_step_diff_eq( nabla_h, Vec3::unit()*EPS0(), // dE/dt coeff - sigma + gamma_sum*EPS0(), // E coeff - gamma_sum.elem_mul(sigma) + gamma_prod*EPS0(), // \int E coeff - gamma_prod.elem_mul(sigma), // \int \int E coeff + sigma, // E coeff + Vec3::zero(), // \int E coeff + Vec3::zero(), // \int \int E coeff delta_t, e_prev, - e_int_prev, - e_int_int_prev, + Vec3::zero(), // e_int_prev + Vec3::zero(), // e_int_int_prev ); - // compute coefficients in: - // nabla_h = de_dt_coeff * dE/dt + prev_e_coeff * prev_e + prev_int_e_coeff * prev_int_e + prev_int_int_e_coeff * prev_int_int_e - // let de_dt_coeff = de_dt_coeff + (e_coeff + (int_e_coeff + int_int_e_coeff*THIRD()*delta_t)*HALF()*delta_t)*delta_t; - // let prev_e_coeff = e_coeff + (int_e_coeff + int_int_e_coeff*HALF()*delta_t)*delta_t; - // let prev_int_e_coeff = int_e_coeff + int_int_e_coeff*delta_t; - // let prev_int_int_e_coeff = int_int_e_coeff; - - // // compute coefficients in: - // // nabla_h = de_dt_coeff * dE/dt + offset - // let offset = e_prev.elem_mul(prev_e_coeff) - // + e_int_prev.elem_mul(prev_int_e_coeff) - // + e_int_int_prev.elem_mul(prev_int_int_e_coeff); - - // let de_dt = (nabla_h - offset).elem_div(de_dt_coeff); - Cell { state: CellState { e: e_prev + de_dt*twice_delta_t, - e_integral: e_int_prev + (e_prev + de_dt*delta_t)*twice_delta_t, - e_integral_integral: e_int_int_prev + (e_int_prev + (e_prev + de_dt*THIRD()*twice_delta_t)*HALF()*twice_delta_t)*twice_delta_t, ..self.state }, mat: self.mat, @@ -908,12 +808,6 @@ impl Cell { pub struct CellState { e: Vec3, h: Vec3, - // The following are used to achieve stretched coordinates (useful for PML). - // XXX: not all materials need these: moving it onto the material could reduce memory footprint - b_integral: Vec3, - h_integral: Vec3, - e_integral: Vec3, - e_integral_integral: Vec3, } impl CellState { @@ -1277,7 +1171,7 @@ mod test { let b = boundary_ness.elem_pow(3.0); let coord_stretch = b * (0.5 / timestep); Static { - coord_stretch, ..Default::default() + conductivity: coord_stretch, ..Default::default() // TODO PML coord_stretch } }); state @@ -1440,7 +1334,7 @@ mod test { _ => unreachable!(), }.into(); Static { - coord_stretch, ..Default::default() + conductivity: coord_stretch, ..Default::default() // TODO PML coord_stretch } }); state @@ -1458,7 +1352,7 @@ mod test { // permute the stretching so that it shouldn't interfere with the wave let coord_stretch = shuffle(cs); Static { - coord_stretch, ..Default::default() + conductivity: coord_stretch, ..Default::default() // TODO PML coord_stretch } }); let center = Index(state.size().0 / 2);