diff --git a/src/mat.rs b/src/mat.rs index 80d2c53..2aa80f9 100644 --- a/src/mat.rs +++ b/src/mat.rs @@ -1,6 +1,7 @@ use crate::{CellState, consts}; use crate::flt::{Flt, Real}; use crate::geom::{Line2d, Vec2, Vec3, Polygon2d}; +use crate::sim::{UpdateMethod, UpdateMethodMut}; use lazy_static::lazy_static; use log::{debug, trace}; @@ -10,10 +11,9 @@ use std::cmp::Ordering; #[enum_dispatch] pub trait Material { - /// Return \sigma, the electrical conductivity. - /// For a vacuum, this is zero. For a perfect conductor, \inf. - fn conductivity(&self) -> Vec3 { - Default::default() + fn update_method_mut<'a>(&'a mut self) -> UpdateMethodMut<'a> { + // by default, behave as a vacuum + UpdateMethodMut::Maxwell { conductivity: Vec3::zero() } } /// Return the magnetization. fn m(&self) -> Vec3 { @@ -27,12 +27,27 @@ pub trait Material { } pub trait MaterialExt { + fn update_method<'a>(&'a self) -> UpdateMethod<'a>; + fn conductivity(&self) -> Vec3; fn is_vacuum(&self) -> bool; } impl MaterialExt for M { + fn update_method<'a>(&'a self) -> UpdateMethod<'a> { + match unsafe { &mut *(self as *const M as *mut M) }.update_method_mut() { + UpdateMethodMut::Maxwell { conductivity } => UpdateMethod::Maxwell { conductivity}, + UpdateMethodMut::Pml { pml, pseudo_conductivity } => + UpdateMethod::Pml { pml, pseudo_conductivity }, + } + } + fn conductivity(&self) -> Vec3 { + match self.update_method() { + UpdateMethod::Maxwell { conductivity } => conductivity, + UpdateMethod::Pml { pseudo_conductivity, .. } => pseudo_conductivity, + } + } fn is_vacuum(&self) -> bool { - self.conductivity() == Default::default() + self.conductivity() == Vec3::zero() } } @@ -54,8 +69,8 @@ impl Static { } impl Material for Static { - fn conductivity(&self) -> Vec3 { - self.conductivity + fn update_method_mut<'a>(&'a mut self) -> UpdateMethodMut<'a> { + UpdateMethodMut::Maxwell { conductivity: self.conductivity } } fn m(&self) -> Vec3 { self.m @@ -89,8 +104,8 @@ impl Conductor { } impl Material for Conductor { - fn conductivity(&self) -> Vec3 { - self.conductivity + fn update_method_mut<'a>(&'a mut self) -> UpdateMethodMut<'a> { + UpdateMethodMut::Maxwell { conductivity: self.conductivity } } } @@ -102,9 +117,10 @@ pub trait PiecewiseLinearFerromagnet { } impl Material for T { - fn conductivity(&self) -> Vec3 { + fn update_method_mut<'a>(&'a mut self) -> UpdateMethodMut<'a> { let c = T::conductivity(); - Vec3::new(c, c, c) + let conductivity = Vec3::new(c, c, c); + UpdateMethodMut::Maxwell { conductivity } } fn m(&self) -> Vec3 { self.m() diff --git a/src/sim.rs b/src/sim.rs index 8947e60..1037590 100644 --- a/src/sim.rs +++ b/src/sim.rs @@ -1,6 +1,6 @@ use crate::{flt::{Flt, Real}, consts}; use crate::geom::{Coord, Cube, Index, InvertedRegion, Meters, Region, Vec3, Vec3u}; -use crate::mat::{self, GenericMaterial, Material}; +use crate::mat::{self, GenericMaterial, Material, MaterialExt as _}; use crate::stim::AbstractStimulus; use dyn_clone::{self, DynClone}; use log::trace; @@ -13,6 +13,43 @@ use std::iter::Sum; pub type StaticSim = SimState; +#[derive(Default, Copy, Clone, PartialEq)] +pub struct PmlState { + /// Auxiliary field used when stepping E + aux_e: Vec3, + /// Auxiliary field used when stepping B + aux_b: Vec3, +} + +impl PmlState { + pub fn new() -> Self { + Self::default() + } +} + +/// Used to instruct the Sim how to advance the state at a given location. +/// +/// `UpdateMethod::Maxwell` means to apply Maxwell's equations as usual, with +/// the provided parameters. +/// +/// `UpdateMethod::Pml` means to apply a modified formulation of Maxwell's equations +/// which are reflectionless for normally incident waves. This allows the material to behave +/// as a "Perfectly Matched Layer" (PML), dissipating energy AS IF it was a conductor +/// with the provided `pseudo_conductivity`. +#[derive(PartialEq)] +pub enum UpdateMethodMut<'a> { + /// `conductivity` is traditionally represented by $\sigma$ + Maxwell { conductivity: Vec3 }, + Pml { pml: &'a mut PmlState, pseudo_conductivity: Vec3 }, +} + +#[derive(PartialEq)] +pub enum UpdateMethod<'a> { + /// `conductivity` is traditionally represented by $\sigma$ + Maxwell { conductivity: Vec3 }, + Pml { pml: &'a PmlState, pseudo_conductivity: Vec3 }, +} + pub trait GenericSim: Send + Sync + DynClone { fn sample(&self, pos: Meters) -> Cell; fn impulse_e_meters(&mut self, pos: Meters, amount: Vec3); @@ -580,6 +617,7 @@ impl Cell { } pub fn current_density(&self) -> Vec3 { + // TODO: does this make sense for Pml? let conductivity = self.mat.conductivity(); self.e().elem_mul(conductivity) } @@ -664,7 +702,6 @@ impl Cell { let twice_delta_t = TWO() * delta_t; let half_delta_t = HALF() * delta_t; - let h_prev = self.h(); let b_prev = self.b(); let (delta_ey_x, delta_ez_x) = match (neighbors.right, neighbors.left) { @@ -691,28 +728,33 @@ impl Cell { // \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 - Vec3::zero(), // B coeff - Vec3::zero(), // \int B coeff - Vec3::zero(), // \int \int B coeff - delta_t, - b_prev, - Vec3::zero(), // b_int_prev, - Vec3::zero(), // b_int_int_prev - ); + match self.mat.update_method() { + UpdateMethod::Maxwell { .. } => { + // From: $\nabla x E = -dB/dt$ + let db_dt = solve_step_diff_eq( + 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, + Vec3::zero(), // b_int_prev, + Vec3::zero(), // b_int_int_prev + ); - // Maxwell's equation gave us delta_b. Material parameters let us turn that into h. - let delta_b = db_dt * twice_delta_t; - let h = self.mat.step_h(&self.state, delta_b); - Cell { - state: CellState { - h: h.into(), - ..self.state - }, - mat: self.mat, + // Maxwell's equation gave us delta_b. Material parameters let us turn that into h. + let delta_b = db_dt * twice_delta_t; + let h = self.mat.step_h(&self.state, delta_b); + Cell { + state: CellState { + h: h.into(), + ..self.state + }, + mat: self.mat, + } + } + UpdateMethod::Pml { .. } => unimplemented!(), } } @@ -751,9 +793,6 @@ impl Cell { use consts::real::{EPS0, HALF, THIRD, TWO, ZERO}; let twice_delta_t = TWO() * delta_t; - let sigma = self.mat.conductivity(); - - let e_prev = self.e(); let (delta_hy_x, delta_hz_x) = match (neighbors.left, neighbors.right) { @@ -780,26 +819,32 @@ impl Cell { // \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)); + match self.mat.update_method() { + UpdateMethod::Maxwell { conductivity } => { + let sigma = conductivity; - // 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, // E coeff - Vec3::zero(), // \int E coeff - Vec3::zero(), // \int \int E coeff - delta_t, - e_prev, - Vec3::zero(), // e_int_prev - Vec3::zero(), // e_int_int_prev - ); + // 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, // E coeff + Vec3::zero(), // \int E coeff + Vec3::zero(), // \int \int E coeff + delta_t, + e_prev, + Vec3::zero(), // e_int_prev + Vec3::zero(), // e_int_int_prev + ); - Cell { - state: CellState { - e: e_prev + de_dt*twice_delta_t, - ..self.state - }, - mat: self.mat, + Cell { + state: CellState { + e: e_prev + de_dt*twice_delta_t, + ..self.state + }, + mat: self.mat, + } + } + UpdateMethod::Pml { .. } => unimplemented!(), } } }