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

@@ -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()
         }

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<M: Material> 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<T> From<T> for Static

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;

src/sim.rs

@@ -516,59 +516,6 @@ impl<M> Cell<M> {
         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<M: Material> Cell<M> {
         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.ey() - self.ey()),
-                inv_feature_size * ((right.ez() - self.ez()) + gamma.z() * (right.int_ez() - self.int_ez())),
+                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.ex() - self.ex()),
-                inv_feature_size * ((down.ez() - self.ez()) + gamma.z() * (down.int_ez() - self.int_ez())),
+                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_.ex() - self.ex()),
-                inv_feature_size * ((in_.ey() - self.ey()) + gamma.y() * (in_.int_ey() - self.int_ey())),
+                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,
+            nabla_e,
-            Vec3::unit(),  // dB/dt coeff
+            -Vec3::unit(),  // dB/dt coeff
-            gamma_sum,    // B coeff
+            Vec3::zero(),    // B coeff
-            gamma_prod,   // \int 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<M: Material> Cell<M> {
         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<M: Material> Cell<M> {
 
         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.hy() - left.hy()),
-                inv_feature_size * ((self.hz() - left.hz()) + gamma.z() * (self.int_hz() - left.int_hz())),
+                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.hx() - up.hx()),
-                inv_feature_size * ((self.hz() - up.hz()) + gamma.z() * (self.int_hz() - up.int_hz())),
+                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.hx() - out.hx()),
-                inv_feature_size * ((self.hy() - out.hy()) + gamma.y() * (self.int_hy() - out.int_hy())),
+                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
+            sigma, // E coeff
-            gamma_sum.elem_mul(sigma) + gamma_prod*EPS0(), // \int E coeff
+            Vec3::zero(), // \int E coeff
-            gamma_prod.elem_mul(sigma), // \int \int E coeff
+            Vec3::zero(), // \int \int E coeff
             delta_t,
             e_prev,
-            e_int_prev,
+            Vec3::zero(), // e_int_prev
-            e_int_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<M: Material> Cell<M> {
 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);