From f6e4ed8cb9666e70094aef9f94a74c337ca1266b Mon Sep 17 00:00:00 2001
From: Colin <colin@uninsane.org>
Date: Sat, 19 Sep 2020 20:02:11 -0700
Subject: [PATCH] Add a UPML boundary condition (needs improving in the
 corners, though)

---
 examples/boundary_conditions.rs | 25 +++++++++++++++++++++++++
 src/driver.rs                   | 24 ++++++++++++++++++++++++
 src/mat.rs                      |  5 ++++-
 src/render.rs                   |  2 +-
 src/sim.rs                      |  8 ++++----
 5 files changed, 58 insertions(+), 6 deletions(-)
 create mode 100644 examples/boundary_conditions.rs

diff --git a/examples/boundary_conditions.rs b/examples/boundary_conditions.rs
new file mode 100644
index 0000000..9c463fa
--- /dev/null
+++ b/examples/boundary_conditions.rs
@@ -0,0 +1,25 @@
+use coremem::{Driver, mat};
+
+fn main() {
+    coremem::init_logging();
+    let width = 201;
+    let boundary = 20;
+    let mut driver = Driver::new(width, width, 1e-3 /* feature size */);
+    driver.add_y4m_renderer(&*format!("boundary_conditions_upml-{}.y4m", boundary));
+    // driver.add_term_renderer();
+
+    //driver.add_boundary(boundary, 0.1);
+    driver.add_upml_boundary(boundary);
+
+    loop {
+        let imp = if driver.state.step_no() < 50 { 
+            30000.0 * ((driver.state.step_no() as f64)*0.04*std::f64::consts::PI).sin()
+        } else {
+            0.0
+        };
+        for y in 0..width {
+            driver.state.impulse_bz(100, y, (imp / 3.0e8) as _);
+        }
+        driver.step();
+    }
+}
diff --git a/src/driver.rs b/src/driver.rs
index 8c93ff6..58e1132 100644
--- a/src/driver.rs
+++ b/src/driver.rs
@@ -1,7 +1,9 @@
 use crate::{flt::Flt, mat};
+use crate::consts;
 use crate::meas::{self, AbstractMeasurement};
 use crate::render::{self, MultiRenderer, Renderer};
 use crate::sim::SimState;
+use crate::vec3::Vec3;
 
 use log::{info, trace};
 use std::path::PathBuf;
@@ -72,6 +74,28 @@ impl Driver {
         }
     }
 
+    pub fn add_upml_boundary(&mut self, thickness: u32) {
+        // Based on explanation here (slide 63): https://empossible.net/wp-content/uploads/2020/01/Lecture-The-Perfectly-Matched-Layer.pdf
+        for inset in 0..thickness {
+            let depth = thickness - inset;
+            let conductivity = 0.5 * consts::EPS0 / self.state.timestep() * (depth as Flt/thickness as Flt).powf(3.0);
+            for x in inset..self.state.width() - inset {
+                let conductor = mat::Static::anisotropic_conductor(Vec3::new(0.0, conductivity, 0.0));
+                // left
+                *self.state.get_mut(x, inset).mat_mut() = conductor.clone().into();
+                // right
+                *self.state.get_mut(x, self.state.height() - 1 - inset).mat_mut() = conductor.into();
+            }
+            for y in inset..self.state.height() - inset {
+                let conductor = mat::Static::anisotropic_conductor(Vec3::new(conductivity, 0.0, 0.0));
+                // top
+                *self.state.get_mut(inset, y).mat_mut() = conductor.clone().into();
+                // bottom
+                *self.state.get_mut(self.state.width() - 1 - inset, y).mat_mut() = conductor.into();
+            }
+        }
+    }
+
     pub fn step(&mut self) {
         if self.state.step_no() % self.steps_per_frame == 0 {
             trace!("render begin");
diff --git a/src/mat.rs b/src/mat.rs
index dddede7..ab33ef3 100644
--- a/src/mat.rs
+++ b/src/mat.rs
@@ -32,8 +32,11 @@ pub struct Static {
 
 impl Static {
     pub fn conductor(conductivity: Flt) -> Self {
+        Self::anisotropic_conductor(Vec3::unit() * conductivity)
+    }
+    pub fn anisotropic_conductor(conductivity: Vec3) -> Self {
         Self {
-            conductivity: Vec3::unit() * conductivity,
+            conductivity,
             ..Default::default()
         }
     }
diff --git a/src/render.rs b/src/render.rs
index 7542c23..43486ce 100644
--- a/src/render.rs
+++ b/src/render.rs
@@ -63,7 +63,7 @@ impl<'a> RenderSteps<'a> {
         let mut me = Self::new(state, measurements);
         me.render_scalar_field(10.0,  false, 2, |cell| cell.mat().conductivity().mag());
         me.render_scalar_field(100.0, true,  0, |cell| cell.mat().m().mag());
-        if false {
+        if true {
             me.render_b_z_field();
             me.render_e_xy_field();
         } else {
diff --git a/src/sim.rs b/src/sim.rs
index 878e2f8..c5f0d97 100644
--- a/src/sim.rs
+++ b/src/sim.rs
@@ -32,7 +32,7 @@ impl<M: Material + Default> SimState<M> {
 impl<M: Material + Clone + Default + Send + Sync> SimState<M> {
     pub fn step(&mut self) {
         use consts::real::*;
-        let half_time_step = HALF() * self.timestep();
+        let half_time_step = HALF() * Real::from_inner(self.timestep());
         let w = self.width() as usize;
         let h = self.height() as usize;
         let feature_size = self.feature_size;
@@ -113,7 +113,7 @@ impl<M: Material> SimState<M> {
 
 impl<M> SimState<M> {
     pub fn time(&self) -> Flt {
-        (self.timestep() * self.step_no as Flt).into()
+        self.timestep() * self.step_no as Flt
     }
 
     pub fn step_no(&self) -> u64 {
@@ -147,8 +147,8 @@ impl<M> SimState<M> {
         &mut self.cells[[y as _, x as _]]
     }
 
-    fn timestep(&self) -> Real {
-        self.feature_size / consts::real::C()
+    pub fn timestep(&self) -> Flt {
+        (self.feature_size / consts::real::C()).into()
     }
 }