colin/fdtd-coremem
1
0
Fork 0
Code Issues Pull Requests Projects Releases Wiki Activity

Fix some issues related to the Courant number

Browse Source 
step_h was using the wrong timestep; fixing it caused update issues
based on the Courant number being 1.0, i.e. too large. Reducing it to
\sqrt{3}/3 fixes stability issues. Unfortunately, a non-Unity Courant
value prevents me from using the ideal Absorbing Boundary Condition I
had just implemented, so that's been disabled.
This commit is contained in:
Colin
2020-10-03 19:39:45 -07:00
parent c33317c805
commit 0642f52172
4 changed files with 158 additions and 63 deletions
Download Patch File Download Diff File Expand all files Collapse all files

22
examples/em_reflection.rs
View File

@@ -3,13 +3,14 @@ use coremem::geom::Index;
fn main() { fn main() {
coremem::init_logging(); coremem::init_logging();
let width = 201; let width = 401;
let size = Index((width, 201, 1).into()); let height = 401;
let size = Index((width, height, 1).into());
let mut driver = Driver::new(size, 1e-6 /* feature size */); let mut driver = Driver::new(size, 1e-6 /* feature size */);
driver.add_y4m_renderer("em_reflection.y4m"); driver.add_y4m_renderer("em_reflection.y4m");
// driver.add_term_renderer(); // driver.add_term_renderer();
driver.add_upml_boundary(Index((20, 40, 0).into())); // driver.add_upml_boundary(Index((20, 40, 0).into()));
// for y in 75..100 { // for y in 75..100 {
// for x in 0..width { // for x in 0..width {
// // from https://www.thoughtco.com/table-of-electrical-resistivity-conductivity-608499 // // from https://www.thoughtco.com/table-of-electrical-resistivity-conductivity-608499
@@ -23,18 +24,27 @@ fn main() {
loop { loop {
let imp = if driver.state.step_no() < 50 { let imp = if driver.state.step_no() < 50 {
3e4 * ((driver.state.step_no() as f64)*0.04*std::f64::consts::PI).sin() 1e4 * ((driver.state.step_no() as f64)*0.04*std::f64::consts::PI).sin()
} else { } else {
0.0 0.0
}; };
// let t = 4e12 * driver.state.time();
// let imp = 1e4 * (-t*t).exp();
// state.impulse_ex(50, 50, imp); // state.impulse_ex(50, 50, imp);
// state.impulse_ey(50, 50, imp); // state.impulse_ey(50, 50, imp);
// state.impulse_bz(20, 20, (imp / 3.0e8) as _); // state.impulse_bz(20, 20, (imp / 3.0e8) as _);
// state.impulse_bz(80, 20, (imp / 3.0e8) as _); // state.impulse_bz(80, 20, (imp / 3.0e8) as _);
for x in 0..width { for y in height/4..height*3/4 {
let loc = Index((x, 50, 0).into()); //for y in 0..height {
let loc = Index((200, y, 0).into());
driver.dyn_state().impulse_ez(loc, imp as _); driver.dyn_state().impulse_ez(loc, imp as _);
driver.dyn_state().impulse_hy(loc, (imp/376.730) as _);
} }
// for x in 0..width {
// let loc = Index((x, 200, 0).into());
// driver.dyn_state().impulse_ez(loc, imp as _);
// driver.dyn_state().impulse_hx(loc, (imp/376.730) as _);
// }
driver.step(); driver.step();
//thread::sleep(time::Duration::from_millis(67)); //thread::sleep(time::Duration::from_millis(67));
} }

8
src/lib.rs
View File

@@ -36,10 +36,10 @@ pub(crate) mod flt {
#[cfg(debug)] #[cfg(debug)]
pub(crate) mod defaults { pub(crate) mod defaults {
use super::*; use super::*;
pub(crate) type Real = R32; //pub(crate) type Real = R32;
pub type Flt = f32; //pub type Flt = f32;
//pub(crate) type Real = R64; pub(crate) type Real = R64;
//pub type Flt = f64; pub type Flt = f64;
} }
#[cfg(not(debug))] #[cfg(not(debug))]
pub(crate) mod defaults { pub(crate) mod defaults {

19
src/render.rs
View File

@@ -5,7 +5,7 @@ use crate::mat;
use crate::sim::{Cell, GenericSim}; use crate::sim::{Cell, GenericSim};
use crate::meas::AbstractMeasurement; use crate::meas::AbstractMeasurement;
use font8x8::{BASIC_FONTS, GREEK_FONTS, UnicodeFonts as _}; use font8x8::{BASIC_FONTS, GREEK_FONTS, UnicodeFonts as _};
use log::trace; use log::{trace, info};
use image::{RgbImage, Rgb}; use image::{RgbImage, Rgb};
use imageproc::{pixelops, drawing}; use imageproc::{pixelops, drawing};
use std::fs::File; use std::fs::File;
@@ -62,8 +62,14 @@ struct RenderSteps<'a> {
impl<'a> RenderSteps<'a> { impl<'a> RenderSteps<'a> {
fn render(state: &'a dyn GenericSim, measurements: &'a [Box<dyn AbstractMeasurement>], z: u32) -> RgbImage { fn render(state: &'a dyn GenericSim, measurements: &'a [Box<dyn AbstractMeasurement>], z: u32) -> RgbImage {
let width = 768; let mut width = 1920;
let height = width * state.height() / state.width(); let mut height = width * state.height() / state.width();
if height > 1080 {
let stretch = 1080 as f32 / height as f32;
width = (width as f32 * stretch) as _;
height = 1080;
}
trace!("rendering at {}x{}", width, height);
let mut me = Self::new(state, measurements, width, height, z); let mut me = Self::new(state, measurements, width, height, z);
me.render_scalar_field(10.0, false, 2, |cell| cell.mat().conductivity().mag()); 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()); me.render_scalar_field(100.0, true, 0, |cell| cell.mat().m().mag());
@@ -161,9 +167,10 @@ impl<'a> RenderSteps<'a> {
for x in 0..8 { for x in 0..8 {
if (bmp & 1 << x) != 0 { if (bmp & 1 << x) != 0 {
let real_x = 2 + i as u32*8 + x; let real_x = 2 + i as u32*8 + x;
let real_y = y as u32 + self.im.height() - 10 - meas_no as u32 * 8; if let Some(real_y) = (y as u32 + self.im.height()).checked_sub(10 + meas_no as u32 * 8) {
if real_x < self.im.width() { if real_x < self.im.width() {
self.im.put_pixel(real_x, real_y, Rgb([0, 0, 0])); self.im.put_pixel(real_x, real_y, Rgb([0, 0, 0]));
}
} }
} }
} }

172
src/sim.rs
View File

@@ -11,6 +11,7 @@ use std::iter::Sum;
pub trait GenericSim { pub trait GenericSim {
fn sample(&self, pos: Meters) -> Cell<mat::Static>; fn sample(&self, pos: Meters) -> Cell<mat::Static>;
fn impulse_e_meters(&mut self, pos: Meters, amount: Vec3); fn impulse_e_meters(&mut self, pos: Meters, amount: Vec3);
fn impulse_h_meters(&mut self, pos: Meters, amount: Vec3);
fn size(&self) -> Index; fn size(&self) -> Index;
fn feature_size(&self) -> Flt; fn feature_size(&self) -> Flt;
fn timestep(&self) -> Flt; fn timestep(&self) -> Flt;
@@ -55,6 +56,10 @@ impl<'a> dyn GenericSim + 'a {
self.impulse_e_meters(pos.to_meters(self.feature_size()), amt) self.impulse_e_meters(pos.to_meters(self.feature_size()), amt)
} }
pub fn impulse_h<C: Coord>(&mut self, pos: C, amt: Vec3) {
self.impulse_h_meters(pos.to_meters(self.feature_size()), amt)
}
pub fn impulse_ex<C: Coord>(&mut self, c: C, ex: Flt) { pub fn impulse_ex<C: Coord>(&mut self, c: C, ex: Flt) {
self.impulse_e(c, Vec3::new(ex, 0.0, 0.0)); self.impulse_e(c, Vec3::new(ex, 0.0, 0.0));
} }
@@ -64,6 +69,17 @@ impl<'a> dyn GenericSim + 'a {
pub fn impulse_ez<C: Coord>(&mut self, c: C, ez: Flt) { pub fn impulse_ez<C: Coord>(&mut self, c: C, ez: Flt) {
self.impulse_e(c, Vec3::new(0.0, 0.0, ez)); self.impulse_e(c, Vec3::new(0.0, 0.0, ez));
} }
pub fn impulse_hx<C: Coord>(&mut self, c: C, hx: Flt) {
self.impulse_h(c, Vec3::new(hx, 0.0, 0.0));
}
pub fn impulse_hy<C: Coord>(&mut self, c: C, hy: Flt) {
self.impulse_h(c, Vec3::new(0.0, hy, 0.0));
}
pub fn impulse_hz<C: Coord>(&mut self, c: C, hz: Flt) {
self.impulse_h(c, Vec3::new(0.0, 0.0, hz));
}
} }
#[derive(Default)] #[derive(Default)]
@@ -88,7 +104,8 @@ impl<M: Material + Default> SimState<M> {
impl<M: Material + Clone + Default + Send + Sync> SimState<M> { impl<M: Material + Clone + Default + Send + Sync> SimState<M> {
pub fn step(&mut self) { pub fn step(&mut self) {
use consts::real::*; use consts::real::*;
let half_time_step = HALF() * Real::from_inner(self.timestep()); let time_step = Real::from_inner(self.timestep());
let half_time_step = HALF() * time_step;
let feature_size = self.feature_size; let feature_size = self.feature_size;
let mut scratch_cells = std::mem::replace(&mut self.scratch, Default::default()); let mut scratch_cells = std::mem::replace(&mut self.scratch, Default::default());
@@ -138,6 +155,10 @@ impl<M: Material> GenericSim for SimState<M> {
self.get_mut(pos).state.e += amount; self.get_mut(pos).state.e += amount;
} }
fn impulse_h_meters(&mut self, pos: Meters, amount: Vec3) {
self.get_mut(pos).state.h += amount;
}
fn size(&self) -> Index { fn size(&self) -> Index {
Index(Vec3u::new( Index(Vec3u::new(
self.cells.shape()[2] as _, self.cells.shape()[2] as _,
@@ -149,7 +170,14 @@ impl<M: Material> GenericSim for SimState<M> {
self.feature_size.into() self.feature_size.into()
} }
fn timestep(&self) -> Flt { fn timestep(&self) -> Flt {
(self.feature_size / consts::real::C()).into() // XXX this has to be multiplied by a Courant Number in order to ensure stability.
// For 3d, we need 3 full timesteps in order to communicate information across the corner
// of a grid. The distance traveled during that time is \sqrt{3}. Hence each timestep needs
// to be \sqrt{3}/3, or 0.577.
// In 2d, this would be \sqrt{2}/2.
// It's an upper limit though; it should be safe to go lower.
let courant = 0.577;
(self.feature_size / consts::real::C() * courant).into()
} }
fn step_no(&self) -> u64 { fn step_no(&self) -> u64 {
self.step_no self.step_no
@@ -255,8 +283,8 @@ impl<'a, Cell> Neighbors<'a, Cell> {
let [z, y, x] = location; let [z, y, x] = location;
let left = array.get([z, y, x.wrapping_sub(1)]); let left = array.get([z, y, x.wrapping_sub(1)]);
let right = array.get([z, y, x + 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.wrapping_sub(1), x]);
let up = array.get([z, y + 1, x]); let down = array.get([z, y + 1, x]);
let out = array.get([z.wrapping_sub(1), y, x]); let out = array.get([z.wrapping_sub(1), y, x]);
let in_ = array.get([z + 1, y, x]); let in_ = array.get([z + 1, y, x]);
Self { Self {
@@ -266,30 +294,31 @@ impl<'a, Cell> Neighbors<'a, Cell> {
} }
impl<M> Cell<M> { impl<M> Cell<M> {
pub fn ex(&self) -> Flt {
self.state.e().x()
}
pub fn ey(&self) -> Flt {
self.state.e().y()
}
pub fn ez(&self) -> Flt {
self.state.e().z()
}
pub fn e(&self) -> Vec3 { pub fn e(&self) -> Vec3 {
self.state.e() self.state.e()
} }
pub fn hx(&self) -> Flt { pub fn ex(&self) -> Flt {
self.state.h().x() self.e().x()
} }
pub fn hy(&self) -> Flt { pub fn ey(&self) -> Flt {
self.state.h().y() self.e().y()
} }
pub fn hz(&self) -> Flt { pub fn ez(&self) -> Flt {
self.state.h().z() self.e().z()
} }
pub fn h(&self) -> Vec3 { pub fn h(&self) -> Vec3 {
self.state.h() self.state.h()
} }
pub fn hx(&self) -> Flt {
self.h().x()
}
pub fn hy(&self) -> Flt {
self.h().y()
}
pub fn hz(&self) -> Flt {
self.h().z()
}
pub fn mat(&self) -> &M { pub fn mat(&self) -> &M {
&self.mat &self.mat
} }
@@ -302,6 +331,12 @@ impl<M: Material> Cell<M> {
pub fn b(&self) -> Vec3 { pub fn b(&self) -> Vec3 {
(self.h() + self.mat.m()) * consts::real::MU0() (self.h() + self.mat.m()) * consts::real::MU0()
} }
pub fn bx(&self) -> Flt {
self.b().x()
}
pub fn by(&self) -> Flt {
self.b().y()
}
pub fn bz(&self) -> Flt { pub fn bz(&self) -> Flt {
self.b().z() self.b().z()
} }
@@ -328,33 +363,52 @@ impl<M: Material> Cell<M> {
// //
// For discretization, $dt$ becomes `delta_t` and $dx$, $dy$, $dz$ become `feature_size`. // For discretization, $dt$ becomes `delta_t` and $dx$, $dy$, $dz$ become `feature_size`.
// ``` // ```
use consts::real::ZERO; use consts::real::{ZERO, TWO};
let inv_feature_size = Real::from_inner(1.0) / feature_size; let inv_feature_size = Real::from_inner(1.0) / feature_size;
let (delta_ey_x, delta_ez_x) = match neighbors.right { let (delta_ey_x, delta_ez_x) = match (neighbors.right, neighbors.left) {
None => (ZERO(), ZERO()), (Some(right), _) => (
Some(right) => (
inv_feature_size * (right.ey() - self.ey()), inv_feature_size * (right.ey() - self.ey()),
inv_feature_size * (right.ez() - self.ez()), inv_feature_size * (right.ez() - self.ez()),
) ),
// This assumes a Courant value of unity.
// (None, Some(left)) => (
// // inv_feature_size * (self.ey() - left.ey()),
// // inv_feature_size * (self.ez() - left.ez()),
// Real::from_inner(self.bz() - left.bz()) / (TWO() * delta_t),
// -Real::from_inner(self.by() - left.by()) / (TWO() * delta_t),
// ),
_ => (ZERO(), ZERO()),
}; };
let (delta_ex_y, delta_ez_y) = match neighbors.down { let (delta_ex_y, delta_ez_y) = match (neighbors.down, neighbors.up) {
None => (ZERO(), ZERO()), (Some(down), _) => (
Some(down) => (
inv_feature_size * (down.ex() - self.ex()), inv_feature_size * (down.ex() - self.ex()),
inv_feature_size * (down.ez() - self.ez()), inv_feature_size * (down.ez() - self.ez()),
) ),
// (None, Some(up)) => (
// // inv_feature_size * (self.ex() - up.ex()),
// // inv_feature_size * (self.ez() - up.ez()),
// -Real::from_inner(self.bz() - up.bz()) / (TWO() * delta_t),
// Real::from_inner(self.bx() - up.bx()) / (TWO() * delta_t),
// ),
_ => (ZERO(), ZERO()),
}; };
let (delta_ex_z, delta_ey_z) = match neighbors.in_ { let (delta_ex_z, delta_ey_z) = match (neighbors.in_, neighbors.out) {
None => (ZERO(), ZERO()), (Some(in_), _) => (
Some(in_) => (
inv_feature_size * (in_.ex() - self.ex()), inv_feature_size * (in_.ex() - self.ex()),
inv_feature_size * (in_.ey() - self.ey()), inv_feature_size * (in_.ey() - self.ey()),
) ),
// (None, Some(out)) => (
// // inv_feature_size * (self.ex() - out.ex()),
// // inv_feature_size * (self.ey() - out.ey()),
// Real::from_inner(self.by() - out.by()) / (TWO() * delta_t),
// -Real::from_inner(self.bx() - out.bx()) / (TWO() * delta_t),
// ),
_ => (ZERO(), ZERO()),
}; };
let nabla_e = Vec3::from((delta_ez_y - delta_ey_z, delta_ex_z - delta_ez_x, delta_ey_x - delta_ex_y)); let nabla_e = Vec3::from((delta_ez_y - delta_ey_z, delta_ex_z - delta_ez_x, delta_ey_x - delta_ex_y));
let delta_b = -nabla_e * delta_t; // TODO: should scale delta_t by 2? let delta_b = -nabla_e * TWO() * delta_t; // TODO: should scale delta_t by 2?
let h = self.mat.step_h(&self.state, delta_b); let h = self.mat.step_h(&self.state, delta_b);
Cell { Cell {
state: CellState { state: CellState {
@@ -401,35 +455,59 @@ impl<M: Material> Cell<M> {
let sigma: Vec3 = self.mat.conductivity(); let sigma: Vec3 = self.mat.conductivity();
let inv_feature_size = Real::from_inner(1.0) / feature_size; let inv_feature_size = Real::from_inner(1.0) / feature_size;
let delta_t_eps = delta_t/EPS0();
let (delta_hy_x, delta_hz_x) = match neighbors.left { let (delta_hy_x, delta_hz_x) = match (neighbors.left, neighbors.right) {
None => (ZERO(), ZERO()), (Some(left), _) => (
Some(left) => (
inv_feature_size * (self.hy() - left.hy()), inv_feature_size * (self.hy() - left.hy()),
inv_feature_size * (self.hz() - left.hz()), inv_feature_size * (self.hz() - left.hz()),
) ),
// This assumes a Courant value of unity
// (None, Some(right)) => (
// // inv_feature_size * (right.hy() - self.hy()),
// // inv_feature_size * (right.hz() - self.hz()),
// // At the edge, we want a normal wave to have zero reflection.
// // That means E(0, t+1) = E(1, t), assuming no OTHER excitation.
// // But there could be non-normal waves, so we avoid short-circuiting the Nabla H,
// // and instead set it to a value which would cause E(0, t+1) = E(1, t) should there
// // be no other excitation.
// Real::from_inner(right.ez() - self.ez()) / (TWO()*delta_t_eps),
// -Real::from_inner(right.ey() - self.ey()) / (TWO()*delta_t_eps),
// ),
// (None, None) => (ZERO(), ZERO()),
_ => (ZERO(), ZERO()),
}; };
let (delta_hx_y, delta_hz_y) = match neighbors.up { let (delta_hx_y, delta_hz_y) = match (neighbors.up, neighbors.down) {
None => (ZERO(), ZERO()), (Some(up), _) => (
Some(up) => (
inv_feature_size * (self.hx() - up.hx()), inv_feature_size * (self.hx() - up.hx()),
inv_feature_size * (self.hz() - up.hz()), inv_feature_size * (self.hz() - up.hz()),
) ),
// (None, Some(down)) => (
// // inv_feature_size * (down.hx() - self.hx()),
// // inv_feature_size * (down.hz() - self.hz()),
// -Real::from_inner(down.ez() - self.ez()) / (TWO()*delta_t_eps),
// Real::from_inner(down.ex() - self.ex()) / (TWO()*delta_t_eps),
// ),
_ => (ZERO(), ZERO()),
}; };
let (delta_hx_z, delta_hy_z) = match neighbors.out { let (delta_hx_z, delta_hy_z) = match (neighbors.out, neighbors.in_) {
None => (ZERO(), ZERO()), (Some(out), _) => (
Some(out) => (
inv_feature_size * (self.hx() - out.hx()), inv_feature_size * (self.hx() - out.hx()),
inv_feature_size * (self.hy() - out.hy()), inv_feature_size * (self.hy() - out.hy()),
) ),
// (None, Some(in_)) => (
// // inv_feature_size * (in_.hx() - self.hx()),
// // inv_feature_size * (in_.hy() - self.hy()),
// Real::from_inner(in_.ey() - self.ey()) / (TWO()*delta_t_eps),
// -Real::from_inner(in_.ex() - self.ex()) / (TWO()*delta_t_eps),
// ),
_ => (ZERO(), ZERO()),
}; };
let nabla_h = Vec3::from((delta_hz_y - delta_hy_z, delta_hx_z - delta_hz_x, delta_hy_x - delta_hx_y)); let nabla_h = Vec3::from((delta_hz_y - delta_hy_z, delta_hx_z - delta_hz_x, delta_hy_x - delta_hx_y));
let delta_t_eps = delta_t/EPS0();
let rhs = self.state.e().elem_mul(Vec3::unit() - sigma*delta_t_eps) + nabla_h * (TWO()*delta_t_eps); let rhs = self.state.e().elem_mul(Vec3::unit() - sigma*delta_t_eps) + nabla_h * (TWO()*delta_t_eps);
let e_next = rhs.elem_div(Vec3::unit() + sigma*delta_t_eps); let e_next = rhs.elem_div(Vec3::unit() + sigma*delta_t_eps);
Cell { Cell {
state: CellState { state: CellState {
e: e_next, e: e_next,