colin
8c3a26e798
- this will allow it to be used from within the spirv code. - had to change some coremem code which was previously peering into privates or now-unrestricted constraints. - may need to put the serde stuff behind a feature flag (or force nostd?)
2074 lines
75 KiB
Rust
2074 lines
75 KiB
Rust
use crate::geom::{Coord, Cube, Index, InvertedRegion, Meters, Region, Vec3, Vec3u};
|
|
use crate::mat::{self, GenericMaterial, Material, MaterialExt as _};
|
|
use crate::real::{R32, Real};
|
|
use crate::stim::{AbstractStimulus, NoopStimulus};
|
|
use log::trace;
|
|
use ndarray::{Array3, Zip};
|
|
use rayon::prelude::*;
|
|
use serde::{Serialize, Deserialize};
|
|
use std::convert::From;
|
|
use std::iter::Sum;
|
|
|
|
pub mod spirv;
|
|
pub mod units;
|
|
|
|
pub type StaticSim = SimState<R32, mat::Static<R32>>;
|
|
|
|
#[derive(Default, Copy, Clone, PartialEq, Serialize, Deserialize)]
|
|
pub struct PmlState<R> {
|
|
/// Auxiliary fields used when stepping E
|
|
aux_e: PmlAux<R>,
|
|
/// Auxiliary fields used when stepping H
|
|
aux_h: PmlAux<R>,
|
|
}
|
|
|
|
impl<R: Default> PmlState<R> {
|
|
pub fn new() -> Self {
|
|
Self::default()
|
|
}
|
|
}
|
|
|
|
#[derive(Default, Copy, Clone, PartialEq, Serialize, Deserialize)]
|
|
struct PmlAux<R> {
|
|
/// ```tex
|
|
/// $Sy^-1 \ast dFz/dy$, where F is either E or H.
|
|
/// This is used to advance the X coordinate of the opposite field
|
|
/// ```
|
|
sy_conv_dfz_dy: ConvFlt<R>,
|
|
sz_conv_dfy_dz: ConvFlt<R>,
|
|
/// Used to advance the Y coordinate of the opposite field
|
|
sz_conv_dfx_dz: ConvFlt<R>,
|
|
sx_conv_dfz_dx: ConvFlt<R>,
|
|
/// Used to advance the Z coordinate of the opposite field
|
|
sx_conv_dfy_dx: ConvFlt<R>,
|
|
sy_conv_dfx_dy: ConvFlt<R>,
|
|
}
|
|
|
|
#[derive(Default, Copy, Clone, PartialEq, Serialize, Deserialize)]
|
|
struct ConvFlt<R>(R);
|
|
|
|
impl<R: Real> ConvFlt<R> {
|
|
///```tex
|
|
/// Consider $v(t) = (a \delta(t) + b exp[-c t] u(t)) \ast y(t)$
|
|
/// This method will advance from self = $v(t)$ to self = $v(t + \Delta t)$
|
|
/// ```
|
|
#[allow(unused)]
|
|
fn step(&mut self, a: R, b: R, c: R, y: R, delta_t: R) -> R {
|
|
assert!(c > R::zero(), "ConvFlt::step called without decay: {}", c.to_f64());
|
|
//```tex
|
|
// let $v(t) = a \delta(t) \ast y(t) + m(t)$
|
|
// where $m(t) = b exp[-c t] u(t) \ast y(t)$
|
|
// let $k(t) = b exp[-c t] u(t)$, where "k" is for "Kernel".
|
|
// $m(t) = \int_0^{t} k(T) y(t - T)dT$
|
|
// $m(t) = \int_0^{\Delta t} k(T) y(t - T)dT + \int_{\Delta t}^t k(T) y(t - T)dT$
|
|
// $m(t) = \int_0^{\Delta t} k(T) y(t - T)dT + \int_{\Delta t}^t b exp[-c T] y(t - T)dT$
|
|
// $m(t) = \int_0^{\Delta t} k(T) y(t - T)dT + \int_0^{t-\Delta t} b exp[-c T] exp[-c \Delta t] y(t - T)dT$
|
|
// $m(t) = \int_0^{\Delta t} k(T) y(t - T)dT + exp[-c \Delta t] v_prev$
|
|
//
|
|
// Assume that $y(t)$ is constant for this time-step.
|
|
// Looking just at the integral term:
|
|
// $b y \int_0^{\Delta t} exp[-c T] dT$
|
|
// $= b y [-1/c](exp[-c \Delta t] - 1)$
|
|
// So
|
|
// $m(t) = exp[-c \Delta t] v_prev + -b y/c (exp[-c \Delta t] - 1)$
|
|
//```
|
|
// trace!("pre: {}", self.0);
|
|
let e = (-c * delta_t).exp();
|
|
// assert!(self.0 < 1e50 && self.0 > -1e50, "self={}, a={} b={} c={} y={} dt={} e={}", self.0, a, b, c, y, delta_t, e);
|
|
self.0 = self.0*e - b*y/c * (e - R::one());
|
|
// trace!("post: {}", self.0);
|
|
self.0 + a*y
|
|
}
|
|
|
|
///```tex
|
|
/// Consider $v(t) = (\delta(t) - \sigma exp[-\sigma t] u(t)) \ast y(t)$
|
|
/// This method will advance from self = $v(t)$ to self = $v(t + \Delta t)$
|
|
///
|
|
/// The user should call this with $e = exp[-\sigma \Delta t]$.
|
|
/// This is a more computationally-efficient version of `step` for when $a=1$ and
|
|
/// $b = -c$ (i.e. `dc_shift` = 0).
|
|
/// ```
|
|
fn step_only_sigma(&mut self, e: R, y: R) -> R {
|
|
//```tex
|
|
// From the `step` method above,
|
|
// $v(t) = a \delta(t) \ast y(t) + m(t)$
|
|
// $m(t) = exp[-c \Delta t] v_prev + -b y/c (exp[-c \Delta t] - 1)$
|
|
// where $a = 1$, $c=\sigma$, $b=-\sigma$.
|
|
// So:
|
|
// $m(t) = e v_prev + \sigma/\sigma (e - 1) y$
|
|
// and $v(t) = m(t) + y(t)$
|
|
// Simplify this:
|
|
// $v(t) = e (v_prev + y)$
|
|
// $m(t) = v(t) - y$
|
|
//```
|
|
let r = e*(self.0 + y);
|
|
self.0 = r - y;
|
|
r
|
|
}
|
|
}
|
|
|
|
|
|
/// Used to instruct the Sim how to advance the state at a given location.
|
|
#[derive(Default, PartialEq)]
|
|
pub struct StepParametersMut<'a, R> {
|
|
conductivity: Vec3<R>,
|
|
pml: Option<(&'a mut PmlState<R>, PmlParameters<R>)>,
|
|
}
|
|
|
|
impl<'a, R: Real> StepParametersMut<'a, R> {
|
|
pub fn new<R2: Real>(
|
|
conductivity: Vec3<R2>, pml: Option<(&'a mut PmlState<R>, PmlParameters<R>)>
|
|
) -> Self {
|
|
Self {
|
|
conductivity: conductivity.cast(),
|
|
pml,
|
|
}
|
|
}
|
|
pub fn with_conductivity<R2: Real>(mut self, conductivity: Vec3<R2>) -> Self {
|
|
self.conductivity = conductivity.cast();
|
|
self
|
|
}
|
|
pub fn with_pml(mut self, state: &'a mut PmlState<R>, params: PmlParameters<R>) -> Self {
|
|
self.pml = Some((state, params));
|
|
self
|
|
}
|
|
}
|
|
|
|
// /// Like `PmlParameters`, but with some additional things parameterized.
|
|
// #[derive(Copy, Clone, Default, PartialEq, Serialize, Deserialize)]
|
|
// pub struct CpmlParameters {
|
|
// /// K; this is similar to \epsilon_r. Try setting this to 1.0 at the interface and
|
|
// /// increasing it towards a larger constant at the simulation boundary.
|
|
// propagation_constant: Vec3<Real>,
|
|
// /// sigma; this is the main parameter that should be varied by location within
|
|
// /// the sim space
|
|
// pseudo_conductivity: Vec3<Real>,
|
|
// /// a; this combats artifacting from DC waves. Try setting this to a small constant
|
|
// /// at the interface and decreasing it to zero at the sim boundary.
|
|
// dc_shift: Vec3<Real>,
|
|
// }
|
|
//
|
|
// impl CpmlParameters {
|
|
// pub fn new(pseudo_conductivity: Vec3<Real>, propagation_constant: Vec3<Real>, dc_shift: Vec3<Real>) -> Self {
|
|
// Self {
|
|
// pseudo_conductivity, propagation_constant, dc_shift
|
|
// }
|
|
// }
|
|
// }
|
|
|
|
/// PML applies 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(Copy, Clone, Default, PartialEq, Serialize, Deserialize)]
|
|
pub struct PmlParameters<R> {
|
|
/// sigma; this is the main parameter that should be varied by location within
|
|
/// the sim space. Initialize it to zero at the interface and increase it
|
|
/// to a large number at the simulation boundary.
|
|
pseudo_conductivity: Vec3<R>,
|
|
}
|
|
|
|
impl<R: Real> PmlParameters<R> {
|
|
pub fn new<R2: Real>(pseudo_conductivity: Vec3<R2>) -> Self {
|
|
Self {
|
|
pseudo_conductivity: pseudo_conductivity.cast(),
|
|
}
|
|
}
|
|
}
|
|
|
|
/// See `StepParametersMut`
|
|
#[derive(PartialEq)]
|
|
pub struct StepParameters<'a, R> {
|
|
conductivity: Vec3<R>,
|
|
pml: Option<(&'a PmlState<R>, PmlParameters<R>)>,
|
|
}
|
|
|
|
impl<'a, R: Real> StepParameters<'a, R> {
|
|
pub fn conductivity(&self) -> Vec3<R> {
|
|
self.conductivity
|
|
}
|
|
pub fn pml(&self) -> Option<(&'a PmlState<R>, PmlParameters<R>)> {
|
|
self.pml
|
|
}
|
|
}
|
|
|
|
impl<'a, R> From<StepParametersMut<'a, R>> for StepParameters<'a, R> {
|
|
fn from(p: StepParametersMut<'a, R>) -> Self {
|
|
Self {
|
|
conductivity: p.conductivity,
|
|
pml: p.pml.map(|(s, p)| (&*s, p)),
|
|
}
|
|
}
|
|
}
|
|
|
|
pub trait MaterialSim: GenericSim {
|
|
type Material: PartialEq;
|
|
|
|
fn put_material<C: Coord, M: Into<Self::Material>>(&mut self, pos: C, mat: M);
|
|
// XXX: would ideally return by-ref, but some backends need to return a handle instead
|
|
fn get_material<C: Coord>(&self, pos: C) -> Self::Material;
|
|
|
|
fn fill_region_using<C, Reg, F, M>(&mut self, region: &Reg, f: F)
|
|
where
|
|
Reg: Region,
|
|
F: Fn(C) -> M,
|
|
C: Coord,
|
|
M: Into<Self::Material>
|
|
{
|
|
for z in 0..self.depth() {
|
|
for y in 0..self.height() {
|
|
for x in 0..self.width() {
|
|
let loc = Index((x, y, z).into());
|
|
let meters = loc.to_meters(self.feature_size());
|
|
if region.contains(meters) {
|
|
self.put_material(loc, f(C::from_either(loc, meters)));
|
|
}
|
|
}
|
|
}
|
|
}
|
|
}
|
|
|
|
fn fill_region<Reg: Region, M: Into<Self::Material> + Clone>(&mut self, region: &Reg, mat: M) {
|
|
self.fill_region_using(region, |_idx: Index| mat.clone());
|
|
}
|
|
|
|
fn examine_region<C, Reg, F>(&self, region: &Reg, mut f: F)
|
|
where
|
|
Reg: Region,
|
|
F: FnMut(C, &Self::Material),
|
|
C: Coord
|
|
{
|
|
for z in 0..self.depth() {
|
|
for y in 0..self.height() {
|
|
for x in 0..self.width() {
|
|
let loc = Index((x, y, z).into());
|
|
let meters = loc.to_meters(self.feature_size());
|
|
if region.contains(meters) {
|
|
f(C::from_either(loc, meters), &self.get_material(loc));
|
|
}
|
|
}
|
|
}
|
|
}
|
|
}
|
|
|
|
/// Return true if the given region is filled exclusively with the provided material.
|
|
fn test_region_filled<Reg: Region, M: Into<Self::Material> + Clone>(&self, region: &Reg, mat: M) -> bool {
|
|
let mut all = true;
|
|
self.examine_region(region, |_idx: Index, m: &Self::Material| {
|
|
all = all && m == &mat.clone().into();
|
|
});
|
|
all
|
|
}
|
|
|
|
/// Fill the boundary, where `thickness` describes how far the boundary extends in each
|
|
/// direction, and `f` takes a vec where each coordinate represents how far into the boundary
|
|
/// the location being queried is, in each direction.
|
|
/// e.g. `f((1.0, 0.0, 0.2))` means the location being queried is at either extreme end on the
|
|
/// x axis, is not inside the y axis boundary, and is 20% of the way from the onset of the z
|
|
/// boundary to the edge of the z world.
|
|
fn fill_boundary_using<C, F, M>(&mut self, thickness: C, f: F)
|
|
where
|
|
C: Coord,
|
|
F: Fn(Vec3<f32>) -> M,
|
|
M: Into<Self::Material>,
|
|
{
|
|
// TODO: maybe this function belongs on the Driver?
|
|
let feat = self.feature_size();
|
|
let upper_left = thickness.to_index(feat);
|
|
let size = self.size();
|
|
let lower_right = size - upper_left - Index::new(1, 1, 1);
|
|
let region = InvertedRegion::new(Cube::new(upper_left.to_meters(feat), lower_right.to_meters(feat)));
|
|
self.fill_region_using(®ion, |loc: Index| {
|
|
let depth_x = if loc.x() < upper_left.x() {
|
|
(upper_left.x() - loc.x()) as f32 / upper_left.x() as f32
|
|
} else if loc.x() > lower_right.x() {
|
|
(loc.x() - lower_right.x()) as f32 / upper_left.x() as f32
|
|
} else {
|
|
0.0
|
|
};
|
|
let depth_y = if loc.y() < upper_left.y() {
|
|
(upper_left.y() - loc.y()) as f32 / upper_left.y() as f32
|
|
} else if loc.y() > lower_right.y() {
|
|
(loc.y() - lower_right.y()) as f32 / upper_left.y() as f32
|
|
} else {
|
|
0.0
|
|
};
|
|
let depth_z = if loc.z() < upper_left.z() {
|
|
(upper_left.z() - loc.z()) as f32 / upper_left.z() as f32
|
|
} else if loc.z() > lower_right.z() {
|
|
(loc.z() - lower_right.z()) as f32 / upper_left.z() as f32
|
|
} else {
|
|
0.0
|
|
};
|
|
// println!("{} {}", loc, Vec3::new(depth_x, depth_y, depth_z));
|
|
f(Vec3::new(depth_x, depth_y, depth_z))
|
|
});
|
|
}
|
|
}
|
|
|
|
// XXX the Send/Sync bounds here could be removed with some refactoring
|
|
pub trait GenericSim: SampleableSim {
|
|
fn step(&mut self) {
|
|
// XXX: try not to exercise this path! NoopStimulus is probably a lot of waste.
|
|
self.step_multiple(1, &NoopStimulus);
|
|
}
|
|
fn step_multiple<S: AbstractStimulus>(&mut self, num_steps: u32, s: &S);
|
|
|
|
/// DEPRECATED. Use stimulus instead
|
|
fn impulse_e_meters(&mut self, pos: Meters, amount: Vec3<f32>);
|
|
/// DEPRECATED. Use stimulus instead
|
|
fn impulse_h_meters(&mut self, pos: Meters, amount: Vec3<f32>) {
|
|
self.impulse_b_meters(pos, amount * f32::mu0());
|
|
}
|
|
/// DEPRECATED. Use stimulus instead
|
|
fn impulse_b_meters(&mut self, pos: Meters, amount: Vec3<f32>) {
|
|
self.impulse_h_meters(pos, amount * f32::mu0_inv());
|
|
}
|
|
|
|
fn impulse_e<C: Coord>(&mut self, pos: C, amt: Vec3<f32>) {
|
|
self.impulse_e_meters(pos.to_meters(self.feature_size()), amt)
|
|
}
|
|
|
|
fn impulse_h<C: Coord>(&mut self, pos: C, amt: Vec3<f32>) {
|
|
self.impulse_h_meters(pos.to_meters(self.feature_size()), amt)
|
|
}
|
|
|
|
fn impulse_b<C: Coord>(&mut self, pos: C, amt: Vec3<f32>) {
|
|
self.impulse_b_meters(pos.to_meters(self.feature_size()), amt)
|
|
}
|
|
|
|
fn impulse_ex<C: Coord>(&mut self, c: C, ex: f32) {
|
|
self.impulse_e(c, Vec3::new_x(ex));
|
|
}
|
|
fn impulse_ey<C: Coord>(&mut self, c: C, ey: f32) {
|
|
self.impulse_e(c, Vec3::new_y(ey));
|
|
}
|
|
fn impulse_ez<C: Coord>(&mut self, c: C, ez: f32) {
|
|
self.impulse_e(c, Vec3::new_z(ez));
|
|
}
|
|
|
|
fn impulse_hx<C: Coord>(&mut self, c: C, hx: f32) {
|
|
self.impulse_h(c, Vec3::new_x(hx));
|
|
}
|
|
fn impulse_hy<C: Coord>(&mut self, c: C, hy: f32) {
|
|
self.impulse_h(c, Vec3::new_y(hy));
|
|
}
|
|
fn impulse_hz<C: Coord>(&mut self, c: C, hz: f32) {
|
|
self.impulse_h(c, Vec3::new_z(hz));
|
|
}
|
|
|
|
fn impulse_bx<C: Coord>(&mut self, c: C, bx: f32) {
|
|
self.impulse_b(c, Vec3::new_x(bx));
|
|
}
|
|
fn impulse_by<C: Coord>(&mut self, c: C, by: f32) {
|
|
self.impulse_b(c, Vec3::new_y(by));
|
|
}
|
|
fn impulse_bz<C: Coord>(&mut self, c: C, bz: f32) {
|
|
self.impulse_b(c, Vec3::new_z(bz));
|
|
}
|
|
}
|
|
|
|
impl<'a> dyn SampleableSim + 'a {
|
|
pub fn get<C: Coord>(&self, at: C) -> Sample {
|
|
self.sample(at.to_meters(self.feature_size()))
|
|
}
|
|
|
|
/// Apply `F` to each Cell, and sum the results.
|
|
pub fn map_sum<F, Ret>(&self, f: F) -> Ret
|
|
where
|
|
F: Fn(&Sample) -> Ret + Sync,
|
|
Ret: Sum<Ret> + Send,
|
|
{
|
|
self.map_sum_enumerated(|_at: Index, cell| f(cell))
|
|
}
|
|
|
|
pub fn map_sum_enumerated<C, F, Ret>(&self, f: F) -> Ret
|
|
where C: Coord,
|
|
F: Fn(C, &Sample) -> Ret + Sync,
|
|
Ret: Sum<Ret> + Send,
|
|
{
|
|
let (w, h, d) = (self.width(), self.height(), self.depth());
|
|
(0..d).into_par_iter().map(
|
|
|z| (0..h).into_par_iter().map_with(z,
|
|
|&mut z, y| (0..w).into_par_iter().map_with((z, y),
|
|
|&mut (z, y), x|
|
|
{
|
|
let at = Index(Vec3u::new(x, y, z));
|
|
f(C::from_index(at, self.feature_size()), &self.get(at))
|
|
}))).flatten().flatten().sum()
|
|
}
|
|
|
|
pub fn volume_of_region<R: Region + ?Sized>(&self, region: &R) -> u32 {
|
|
self.map_sum_over(region, |_| 1)
|
|
}
|
|
|
|
/// Apply `F` to each Cell, and sum the results.
|
|
pub fn map_sum_over<F, Ret, Reg>(&self, region: &Reg, f: F) -> Ret
|
|
where
|
|
F: Fn(&Sample) -> Ret + Sync,
|
|
Ret: Sum<Ret> + Default + Send,
|
|
Reg: Region + ?Sized
|
|
{
|
|
self.map_sum_over_enumerated(region, |_at: Index, cell| f(cell))
|
|
}
|
|
|
|
pub fn map_sum_over_enumerated<C, F, Ret, Reg>(&self, region: &Reg, f: F) -> Ret
|
|
where C: Coord,
|
|
F: Fn(C, &Sample) -> Ret + Sync,
|
|
Ret: Sum<Ret> + Default + Send,
|
|
Reg: Region + ?Sized,
|
|
{
|
|
let (w, h, d) = (self.width(), self.height(), self.depth());
|
|
(0..d).into_par_iter().map(
|
|
|z| (0..h).into_par_iter().map_with(z,
|
|
|&mut z, y| (0..w).into_par_iter().map_with((z, y),
|
|
|&mut (z, y), x|
|
|
{
|
|
let at = Index(Vec3u::new(x, y, z));
|
|
let meters = at.to_meters(self.feature_size());
|
|
if region.contains(meters) {
|
|
f(C::from_index(at, self.feature_size()), &self.get(at))
|
|
} else {
|
|
Default::default()
|
|
}
|
|
}))).flatten().flatten().sum()
|
|
}
|
|
|
|
pub fn current<C: Coord>(&self, c: C) -> Vec3<f32> {
|
|
self.get(c).current_density() * self.feature_size() * self.feature_size()
|
|
}
|
|
}
|
|
|
|
// XXX the Send/Sync bounds here could be removed with some refactoring
|
|
pub trait SampleableSim: Send + Sync {
|
|
fn sample(&self, pos: Meters) -> Sample;
|
|
|
|
fn size(&self) -> Index;
|
|
fn feature_size(&self) -> f32;
|
|
fn feature_volume(&self) -> f32 {
|
|
let f = self.feature_size();
|
|
f*f*f
|
|
}
|
|
fn volume(&self) -> f32 {
|
|
let s = self.size().to_meters(self.feature_size());
|
|
s.x() * s.y() * s.z()
|
|
}
|
|
fn timestep(&self) -> f32;
|
|
|
|
fn step_no(&self) -> u64;
|
|
|
|
fn width(&self) -> u32 {
|
|
self.size().x()
|
|
}
|
|
fn height(&self) -> u32 {
|
|
self.size().y()
|
|
}
|
|
fn depth(&self) -> u32 {
|
|
self.size().z()
|
|
}
|
|
fn time(&self) -> f32 {
|
|
self.timestep() * self.step_no() as f32
|
|
}
|
|
|
|
/// Take a "snapshot" of the simulation, dropping all material-specific information.
|
|
fn to_static(&self) -> StaticSim {
|
|
let mut state = SimState::new(self.size(), self.feature_size());
|
|
state.step_no = self.step_no();
|
|
Zip::from(ndarray::indices_of(&state.cells)).and(&mut state.e).and(&mut state.h).par_map_assign_into(
|
|
&mut state.cells,
|
|
|(z, y, x), e, h| {
|
|
let idx = Index((x as u32, y as u32, z as u32).into());
|
|
let cell = self.sample(idx.to_meters(self.feature_size()));
|
|
*e = cell.e().cast();
|
|
*h = cell.h().cast();
|
|
mat::Static {
|
|
conductivity: cell.conductivity().cast(),
|
|
m: cell.m().cast(),
|
|
}
|
|
});
|
|
state
|
|
}
|
|
}
|
|
|
|
#[derive(Default, Clone, Serialize, Deserialize)]
|
|
pub struct SimState<R=R32, M=GenericMaterial<R>> {
|
|
cells: Array3<M>,
|
|
e: Array3<Vec3<R>>,
|
|
h: Array3<Vec3<R>>,
|
|
feature_size: R,
|
|
step_no: u64,
|
|
timestep: R,
|
|
}
|
|
|
|
impl<R, M> SimState<R, M> {
|
|
pub fn size(&self) -> Index {
|
|
Index(Vec3u::new(
|
|
self.cells.shape()[2] as _,
|
|
self.cells.shape()[1] as _,
|
|
self.cells.shape()[0] as _,
|
|
))
|
|
}
|
|
pub fn width(&self) -> u32 {
|
|
self.size().x()
|
|
}
|
|
pub fn height(&self) -> u32 {
|
|
self.size().y()
|
|
}
|
|
pub fn depth(&self) -> u32 {
|
|
self.size().z()
|
|
}
|
|
}
|
|
|
|
impl<R: Real, M> SimState<R, M> {
|
|
pub fn feature_size(&self) -> f32 {
|
|
self.feature_size.to_f32()
|
|
}
|
|
pub fn inv_feature_size(&self) -> f32 {
|
|
(R::one() / self.feature_size).to_f32()
|
|
}
|
|
pub fn timestep(&self) -> f32 {
|
|
self.timestep.to_f32()
|
|
}
|
|
}
|
|
|
|
impl<R: Real, M: Default> SimState<R, M> {
|
|
pub fn new(size: Index, feature_size: f32) -> Self {
|
|
// 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;
|
|
let feature_size = R::from_primitive(feature_size);
|
|
let timestep = feature_size / R::c() * R::from_primitive(courant);
|
|
Self {
|
|
cells: Array3::default((size.z() as _, size.y() as _, size.x() as _)),
|
|
e: Array3::default((size.z() as _, size.y() as _, size.x() as _)),
|
|
h: Array3::default((size.z() as _, size.y() as _, size.x() as _)),
|
|
feature_size,
|
|
timestep,
|
|
..Default::default()
|
|
}
|
|
}
|
|
}
|
|
|
|
impl<R: Real, M: Material<R> + Send + Sync> SimState<R, M> {
|
|
pub fn step(&mut self) {
|
|
let time_step = self.timestep;
|
|
let half_time_step = R::half() * time_step;
|
|
let feature_size = self.feature_size;
|
|
let inv_feature_size = R::one() / feature_size;
|
|
|
|
// advance all the electic fields
|
|
trace!("step_e begin");
|
|
let h_field = &self.h;
|
|
Zip::from(ndarray::indices_of(&self.cells)).and(&mut self.cells).and(&mut self.e).par_for_each(
|
|
|(z, y, x), mat, e| {
|
|
let neighbors = Neighbors::new(h_field, [z, y, x]);
|
|
let h = &h_field[[z, y, x]];
|
|
ECell { mat, e, h }.step(neighbors, half_time_step, inv_feature_size);
|
|
});
|
|
trace!("step_e end");
|
|
|
|
// advance all the magnetic fields
|
|
trace!("step_h begin");
|
|
let e_field = &self.e;
|
|
Zip::from(ndarray::indices_of(&self.cells)).and(&mut self.cells).and(&mut self.h).par_for_each(
|
|
|(z, y, x), mat, h| {
|
|
let neighbors = Neighbors::new(e_field, [z, y, x]);
|
|
let e = &e_field[[z, y, x]];
|
|
HCell { mat, h, e }.step(neighbors, half_time_step, inv_feature_size);
|
|
});
|
|
trace!("step_h end");
|
|
|
|
self.step_no += 1;
|
|
}
|
|
|
|
pub fn apply_stimulus<S: AbstractStimulus + ?Sized>(&mut self, stim: &S) {
|
|
let feature_size = self.feature_size();
|
|
|
|
let t_sec = self.time();
|
|
let timestep = self.timestep;
|
|
trace!("apply_stimulus begin");
|
|
Zip::from(ndarray::indices_of(&self.e)).and(&mut self.e).and(&mut self.h).par_for_each(
|
|
|(z, y, x), e, h| {
|
|
let pos_meters = Index((x as u32, y as u32, z as u32).into()).to_meters(feature_size);
|
|
let (density_e, density_h) = stim.at(t_sec, pos_meters);
|
|
let value_e = density_e.cast::<R>() * timestep.cast::<R>();
|
|
let value_h = density_h.cast::<R>() * timestep.cast::<R>();
|
|
*e += value_e;
|
|
*h += value_h;
|
|
});
|
|
trace!("apply_stimulus end");
|
|
}
|
|
}
|
|
|
|
impl<R: Real, M: Material<R> + Clone + Send + Sync + PartialEq> MaterialSim for SimState<R, M> {
|
|
type Material = M;
|
|
|
|
fn put_material<C: Coord, M2: Into<Self::Material>>(&mut self, pos: C, mat: M2) {
|
|
*self.get_mat_mut(pos) = mat.into();
|
|
}
|
|
|
|
fn get_material<C: Coord>(&self, pos: C) -> M {
|
|
self.get_material(pos).clone()
|
|
}
|
|
}
|
|
|
|
impl<R: Real, M: Material<R> + Send + Sync> SampleableSim for SimState<R, M> {
|
|
fn sample(&self, pos: Meters) -> Sample {
|
|
// TODO: smarter sampling than nearest neighbor?
|
|
let pos_sim = pos.to_index(self.feature_size());
|
|
let idx = [pos_sim.z() as usize, pos_sim.y() as _, pos_sim.x() as _];
|
|
match self.cells.get(idx) {
|
|
Some(mat) => Sample {
|
|
state: CellStateWithM {
|
|
e: self.e[idx].cast(),
|
|
h: self.h[idx].cast(),
|
|
m: mat.m().cast(),
|
|
},
|
|
conductivity: mat.conductivity().cast(),
|
|
},
|
|
None => Default::default(),
|
|
}
|
|
}
|
|
|
|
fn size(&self) -> Index {
|
|
Index(Vec3u::new(
|
|
self.cells.shape()[2] as _,
|
|
self.cells.shape()[1] as _,
|
|
self.cells.shape()[0] as _,
|
|
))
|
|
}
|
|
fn feature_size(&self) -> f32 {
|
|
self.feature_size.to_f32()
|
|
}
|
|
fn timestep(&self) -> f32 {
|
|
self.timestep.to_f32()
|
|
}
|
|
fn step_no(&self) -> u64 {
|
|
self.step_no
|
|
}
|
|
}
|
|
|
|
|
|
impl<R: Real, M: Material<R> + Send + Sync> GenericSim for SimState<R, M> {
|
|
fn step_multiple<S: AbstractStimulus>(&mut self, num_steps: u32, s: &S) {
|
|
for _ in 0..num_steps {
|
|
self.apply_stimulus(s);
|
|
SimState::step(self);
|
|
}
|
|
}
|
|
|
|
fn impulse_e_meters(&mut self, pos: Meters, amount: Vec3<f32>) {
|
|
*self.get_e_mut(pos) += amount.cast();
|
|
}
|
|
|
|
fn impulse_h_meters(&mut self, pos: Meters, amount: Vec3<f32>) {
|
|
self.impulse_b_meters(pos, amount * f32::mu0());
|
|
}
|
|
|
|
fn impulse_b_meters(&mut self, pos: Meters, amount: Vec3<f32>) {
|
|
self.get_hcell(pos).impulse_b(amount.cast());
|
|
}
|
|
}
|
|
|
|
#[allow(unused)]
|
|
impl<R: Real, M> SimState<R, M> {
|
|
fn get_material<C: Coord>(&self, c: C) -> &M {
|
|
let at = c.to_index(self.feature_size.cast());
|
|
&self.cells[at.row_major_idx()]
|
|
}
|
|
fn get_mat_mut<C: Coord>(&mut self, c: C) -> &mut M {
|
|
let at = c.to_index(self.feature_size.cast());
|
|
&mut self.cells[at.row_major_idx()]
|
|
}
|
|
|
|
fn get_e<C: Coord>(&self, c: C) -> &Vec3<R> {
|
|
let at = c.to_index(self.feature_size.cast());
|
|
&self.e[at.row_major_idx()]
|
|
}
|
|
fn get_e_mut<C: Coord>(&mut self, c: C) -> &mut Vec3<R> {
|
|
let at = c.to_index(self.feature_size.cast());
|
|
&mut self.e[at.row_major_idx()]
|
|
}
|
|
|
|
fn get_h<C: Coord>(&self, c: C) -> &Vec3<R> {
|
|
let at = c.to_index(self.feature_size.cast());
|
|
&self.h[at.row_major_idx()]
|
|
}
|
|
fn get_h_mut<C: Coord>(&mut self, c: C) -> &mut Vec3<R> {
|
|
let at = c.to_index(self.feature_size.cast());
|
|
&mut self.h[at.row_major_idx()]
|
|
}
|
|
|
|
fn get_hcell<'a, C: Coord>(&'a mut self, c: C) -> HCell<'a, R, M> {
|
|
let at = c.to_index(self.feature_size.cast());
|
|
let idx = at.row_major_idx();
|
|
HCell {
|
|
e: &self.e[idx],
|
|
h: &mut self.h[idx],
|
|
mat: &mut self.cells[idx],
|
|
}
|
|
}
|
|
}
|
|
|
|
/// Conceptually, one cell looks like this (in 2d):
|
|
///
|
|
/// +-------.-------+
|
|
/// | Ex |
|
|
/// | |
|
|
/// | |
|
|
/// .Ey .Bz .
|
|
/// | |
|
|
/// | |
|
|
/// | |
|
|
/// +-------.-------+
|
|
///
|
|
/// Where the right hand rule indicates that positive Bz is pointing into the page, away from the
|
|
/// reader.
|
|
///
|
|
/// The dot on bottom is Ex of the cell at (x, y+1) and the dot on the right is the Ey of the cell at
|
|
/// (x+1, y). The `+` only indicates the corner of the cell -- nothing of interest is measured at
|
|
/// the pluses.
|
|
///
|
|
/// To adapt this to 3d, we keep the above cross-section, and then add another cross-section at a
|
|
/// half cell depth into the page:
|
|
/// 3d version (there may be multiple solutions):
|
|
///
|
|
/// Ez------.-------+
|
|
/// | By |
|
|
/// | |
|
|
/// | |
|
|
/// .Bx .
|
|
/// | |
|
|
/// | |
|
|
/// | |
|
|
/// +-------.-------+
|
|
///
|
|
/// Altogether, each cell is a cube, with the field vectors lying on the surface of the 1/8th of
|
|
/// the cube closest to the asterisk, at (0, 0, 0)
|
|
///
|
|
/// +------------+------------+
|
|
/// |\ \ \
|
|
/// | \ \ \
|
|
/// | \ \ \
|
|
/// | \ Ez \ By \
|
|
/// | +------------+------------+
|
|
/// | |\ \ \ z
|
|
/// + | \ \ \ \
|
|
/// |\ | \ \ \ \
|
|
/// | \ | \ \ Ex \ \
|
|
/// | \ | *------------+------------+ *--- x
|
|
/// | \|Bx | | | |
|
|
/// | + | | | |
|
|
/// | |\ | | | |
|
|
/// + | \ | | | y
|
|
/// \ | \ | | |
|
|
/// \ | \|Ey |Bz |
|
|
/// \ | +------------+------------+
|
|
/// \| | | |
|
|
/// + | | |
|
|
/// \ | | |
|
|
/// \ | | |
|
|
/// \ | | |
|
|
/// \| | |
|
|
/// +------------+------------+
|
|
///
|
|
#[derive(Clone, Default, Serialize, Deserialize)]
|
|
pub struct Sample {
|
|
state: CellStateWithM<f32>,
|
|
conductivity: Vec3<f32>,
|
|
}
|
|
|
|
struct Neighbors<'a, R> {
|
|
left: Option<&'a Vec3<R>>,
|
|
right: Option<&'a Vec3<R>>,
|
|
up: Option<&'a Vec3<R>>,
|
|
down: Option<&'a Vec3<R>>,
|
|
out: Option<&'a Vec3<R>>,
|
|
in_: Option<&'a Vec3<R>>,
|
|
}
|
|
|
|
impl<'a, R> Neighbors<'a, R> {
|
|
fn new(array: &'a Array3<Vec3<R>>, location: [usize; 3]) -> Self {
|
|
let [z, y, x] = location;
|
|
let left = array.get([z, y, x.wrapping_sub(1)]);
|
|
let right = array.get([z, y, x + 1]);
|
|
let up = array.get([z, y.wrapping_sub(1), x]);
|
|
let down = array.get([z, y + 1, x]);
|
|
let out = array.get([z.wrapping_sub(1), y, x]);
|
|
let in_ = array.get([z + 1, y, x]);
|
|
Self {
|
|
left, right, up, down, out, in_
|
|
}
|
|
}
|
|
}
|
|
|
|
impl Sample {
|
|
pub fn e(&self) -> Vec3<f32> {
|
|
self.state.e()
|
|
}
|
|
pub fn ex(&self) -> f32 {
|
|
self.e().x()
|
|
}
|
|
pub fn ey(&self) -> f32 {
|
|
self.e().y()
|
|
}
|
|
pub fn ez(&self) -> f32 {
|
|
self.e().z()
|
|
}
|
|
|
|
pub fn h(&self) -> Vec3<f32> {
|
|
self.state.h()
|
|
}
|
|
pub fn hx(&self) -> f32 {
|
|
self.h().x()
|
|
}
|
|
pub fn hy(&self) -> f32 {
|
|
self.h().y()
|
|
}
|
|
pub fn hz(&self) -> f32 {
|
|
self.h().z()
|
|
}
|
|
|
|
pub fn b(&self) -> Vec3<f32> {
|
|
self.state.b()
|
|
}
|
|
pub fn bx(&self) -> f32 {
|
|
self.b().x()
|
|
}
|
|
pub fn by(&self) -> f32 {
|
|
self.b().y()
|
|
}
|
|
pub fn bz(&self) -> f32 {
|
|
self.b().z()
|
|
}
|
|
|
|
pub fn m(&self) -> Vec3<f32> {
|
|
self.state.m()
|
|
}
|
|
pub fn conductivity(&self) -> Vec3<f32> {
|
|
self.conductivity
|
|
}
|
|
|
|
pub fn current_density(&self) -> Vec3<f32> {
|
|
// TODO: does this make sense for Pml?
|
|
let conductivity = self.conductivity();
|
|
self.e().elem_mul(conductivity)
|
|
}
|
|
|
|
}
|
|
|
|
/// Solve df/dt in:
|
|
/// nabla_g = df_dt_coeff * df/dt + f_coeff * f + f_int_coeff * \int f + f_int_int_coeff * \int \int f
|
|
#[allow(unused)]
|
|
fn solve_step_diff_eq<R: Real>(
|
|
nabla_g: Vec3<R>,
|
|
df_dt_coeff: Vec3<R>,
|
|
f_coeff: Vec3<R>,
|
|
f_int_coeff: Vec3<R>,
|
|
f_int_int_coeff: Vec3<R>,
|
|
delta_t: R,
|
|
f_prev: Vec3<R>,
|
|
f_int_prev: Vec3<R>,
|
|
f_int_int_prev: Vec3<R>,
|
|
) -> Vec3<R> {
|
|
// ```tex
|
|
// $f(t)$ from $(T-\Delta t)$ to $(T+\Delta t)$ is a line. As such:
|
|
// $f(t') = f_prev + t'*db/dt$ where $t'$ is the time since the last sample of $f$, i.e. $f_prev$
|
|
// and
|
|
// $\int f|t' = f_int_prev + t'*f_prev + 1/2 t'^2 df/dt$
|
|
// $\int \int f|t' = f_int_int_prev + t'*f_int_prev + 1/2 t'^2*f_prev + 1/6 t'^3*df/dt$
|
|
// ```
|
|
// Solve coeffs in:
|
|
// nabla_g = df_dt_coeff * df/dt + f_prev_coeff * f_prev
|
|
// + f_int_prev_coeff * f_int_prev + f_int_int_prev_coeff * f_int_int_prev
|
|
let f_int_int_prev_coeff = f_int_int_coeff;
|
|
let f_int_prev_coeff = f_int_coeff + f_int_int_coeff*delta_t;
|
|
let f_prev_coeff = f_coeff + (f_int_coeff +
|
|
f_int_int_coeff*R::half()
|
|
)*delta_t;
|
|
let df_dt_coeff = df_dt_coeff + (f_coeff +
|
|
(f_int_coeff + f_int_int_coeff*R::third()*delta_t)*R::half()*delta_t
|
|
)*delta_t;
|
|
|
|
// Solve coeffs in nabla_g = df_dt_coeff + offset;
|
|
let offset = f_prev_coeff.elem_mul(f_prev) +
|
|
f_int_prev_coeff.elem_mul(f_int_prev) +
|
|
f_int_int_prev_coeff.elem_mul(f_int_int_prev);
|
|
(nabla_g - offset).elem_div(df_dt_coeff)
|
|
}
|
|
|
|
/// Solve df/dt in:
|
|
/// nabla_g = df_dt_coeff * df/dt + f_coeff * f
|
|
/// Lower-order (more efficient) version of `solve_step_diff_eq`.
|
|
fn solve_step_diff_eq_linear<R: Real>(
|
|
nabla_g: Vec3<R>,
|
|
df_dt_coeff: Vec3<R>,
|
|
f_coeff: Vec3<R>,
|
|
delta_t: R,
|
|
f_prev: Vec3<R>,
|
|
) -> Vec3<R> {
|
|
// Solve coeffs in:
|
|
// nabla_g = df_dt_coeff * df/dt + f_prev_coeff * f_prev
|
|
let df_dt_coeff = df_dt_coeff + f_coeff*delta_t;
|
|
|
|
// Solve coeffs in nabla_g = df_dt_coeff + offset;
|
|
let offset = f_coeff.elem_mul(f_prev);
|
|
(nabla_g - offset).elem_div(df_dt_coeff)
|
|
}
|
|
|
|
struct HCell<'a, R, M> {
|
|
mat: &'a mut M,
|
|
h: &'a mut Vec3<R>,
|
|
e: &'a Vec3<R>,
|
|
}
|
|
|
|
struct ECell<'a, R, M> {
|
|
mat: &'a mut M,
|
|
e: &'a mut Vec3<R>,
|
|
h: &'a Vec3<R>,
|
|
}
|
|
|
|
impl<'a, R: Real, M: Material<R>> HCell<'a, R, M> {
|
|
pub fn b(&self) -> Vec3<R> {
|
|
(*self.h + self.mat.m()) * R::mu0()
|
|
}
|
|
|
|
fn impulse_b(&mut self, delta_b: Vec3<R>) {
|
|
let b_prev = self.b();
|
|
let state = CellState {
|
|
e: *self.e,
|
|
h: *self.h,
|
|
};
|
|
self.mat.step_b(&state, delta_b);
|
|
// mu0 (H + M) = B
|
|
// therefore, H = B mu0^-1 - M
|
|
*self.h = (b_prev + delta_b) * R::mu0_inv() - self.mat.m();
|
|
}
|
|
|
|
fn step(&mut self, neighbors: Neighbors<'_, R>, delta_t: R, inv_feature_size: R) {
|
|
// ```tex
|
|
// Maxwell-Faraday equation: $\nabla x E = -dB/dt$ (1)
|
|
//
|
|
// Note: $\nabla x E$ evaluates to:
|
|
// | i j k | [ dEz/dy - dEy/dz ]
|
|
// | | | |
|
|
// | d/dx d/dy d/dz | = | dEx/dz - dEz/dx |
|
|
// | | | |
|
|
// | Ex Ey Ez | [ dEy/dx - dEx/dy ]
|
|
//
|
|
// To accomodate PML materials, we allow coordinate stretching, specifically
|
|
// in the complex domain.
|
|
//
|
|
// For the rest of this analysis, we'll use $\omega$. This comes from a wave of the form
|
|
// $E(t) = A exp(-j\omega t)$ (2)
|
|
// Maxwell's equations are all linear, so we'll do the derivations considering just a
|
|
// single such wave, knowing that what we derive fully generalizes to any $E(t)$.
|
|
// If we do it right, then $\omega$ will drop out from the end result.
|
|
//
|
|
// Consider: $x' = x + jf(x)$.
|
|
// Then $dx' = (1 + jdf/dx)dx$
|
|
// Let $df/dx = \gamma_x/\omega$. This is chosen to eliminate $\omega$ in the end result.
|
|
// Then, $d/dx = 1/(1 + j\gamma_x/\omega) d/dx'$
|
|
// Likewise, $d/dy = 1/(1 + j\gamma_y/\omega) d/dy'$
|
|
// $d/dz = 1/(1 + j\gamma_z/\omega) d/dz'$
|
|
//
|
|
// Let $S = 1 + j\gamma/\omega$, such that $d/dx = Sx^-1 d/dx'$, etc.
|
|
//
|
|
// Then (1) can be expanded to the following matrix multiplication
|
|
// $ [ 0 -Sz^-1 d/dz' Sy^-1 d/dy' ] $
|
|
// $ | Sz^-1 d/dz' 0 -Sx^-1 d/dx' | E = -dB/dt$
|
|
// $ [ -Sy^-1 d/dy' Sx^-1 d/dx' 0 ]$
|
|
//
|
|
// Normalize by multiplying each side by a diagonal matrix:
|
|
//
|
|
// $ [ 0 -Sy d/dz' Sz d/dy' ] [ SySz ] $
|
|
// $ -| Sx d/dz' 0 -Sz d/dx' | E = diag| SxSz | dB/dt (3)$
|
|
// $ [ -Sx d/dy' Sy d/dx' 0 ] [ SxSy ] $
|
|
//
|
|
// The LHS simplifies to:
|
|
// $ -\nabla x (diag(S) E)$
|
|
// $ = -\nabla x (E + diag(\gamma) j/\omega E)$
|
|
//
|
|
// From (2), we know $B$ is of the form $B = A exp(-j\omega t)$
|
|
// Therefore $\int_0^t B(\tau) d\tau = j/\omega [B(t) - B(0)] = j/\omega B(t)$ (we assume $B(0) = 0$)
|
|
// Thus the LHS becomes:
|
|
// $ = -\nabla x (E + diag(\gamma) \int E)$
|
|
//
|
|
// Expand (3), letting $I = j/\omega$:
|
|
// $ |(1 + \gamma_yI)(1 + \gamma_zI)|$
|
|
// $ -\nabla x (E + diag(\gamma) \int E) = diag[(1 + \gamma_xI)(1 + \gamma_zI)] dB/dt$
|
|
// $ |(1 + \gamma_xI)(1 + \gamma_yI)|$
|
|
//
|
|
// $ |\gamma_y + \gamma_z| |\gamma_y\gamma_z|$
|
|
// $ -\nabla x (E + diag(\gamma) \int E) = dB/dt + diag[\gamma_x + \gamma_z] I dB/dt + diag[\gamma_x\gamma_z] I^2 dB/dt$
|
|
// $ |\gamma_x + \gamma_y| |\gamma_x\gamma_y|$
|
|
//
|
|
// As seen above, $I = j/w$ is equivalent to integration:
|
|
// $ |\gamma_y + \gamma_z| |\gamma_y\gamma_z|$
|
|
// $ -\nabla x (E + diag(\gamma) \int E) = dB/dt + diag[\gamma_x + \gamma_z] B + diag[\gamma_x\gamma_z] \int B$
|
|
// $ |\gamma_x + \gamma_y| |\gamma_x\gamma_y|$
|
|
//
|
|
// Finally, discretization. We're expanding about $t = T$, where we know $B(T-\Delta t)$ and $E(T)$.
|
|
// $B(t)$ from $(T-\Delta t)$ to $(T+\Delta t)$ is a line. As such:
|
|
// $B(t') = B_prev + t'*db/dt$ where $t'$ is the time since the last sample of $B$, i.e. $B_prev$
|
|
// and
|
|
// $\int B|t' = B_int_prev + t'*B_prev + 1/2 t'^2 dB/dt
|
|
//
|
|
// ```
|
|
// let inv_feature_size = Real::from_inner(1.0) / feature_size;
|
|
let twice_delta_t = R::two() * delta_t;
|
|
|
|
let b_prev = self.b();
|
|
|
|
let (delta_ey_x, delta_ez_x) = match (neighbors.right, neighbors.left) {
|
|
(Some(right), _) => (
|
|
inv_feature_size * (right.y() - self.e.y()),
|
|
inv_feature_size * (right.z() - self.e.z()),
|
|
),
|
|
_ => (R::zero(), R::zero()),
|
|
};
|
|
let (delta_ex_y, delta_ez_y) = match (neighbors.down, neighbors.up) {
|
|
(Some(down), _) => (
|
|
inv_feature_size * (down.x() - self.e.x()),
|
|
inv_feature_size * (down.z() - self.e.z()),
|
|
),
|
|
_ => (R::zero(), R::zero()),
|
|
};
|
|
let (delta_ex_z, delta_ey_z) = match (neighbors.in_, neighbors.out) {
|
|
(Some(in_), _) => (
|
|
inv_feature_size * (in_.x() - self.e.x()),
|
|
inv_feature_size * (in_.y() - self.e.y()),
|
|
),
|
|
_ => (R::zero(), R::zero()),
|
|
};
|
|
|
|
let StepParametersMut { conductivity: _, pml } = self.mat.step_parameters_mut();
|
|
// XXX conductivity is unused when stepping h.
|
|
let nabla_e = match pml {
|
|
None => {
|
|
// \nabla x E
|
|
Vec3::from((delta_ez_y - delta_ey_z, delta_ex_z - delta_ez_x, delta_ey_x - delta_ex_y))
|
|
}
|
|
Some((pml, PmlParameters { pseudo_conductivity })) => {
|
|
//```tex
|
|
// S = propagation_const + pseudo_conductivity / (dc_shift + j\omega \mu_0)
|
|
// $S = 1 + sigma/(j\omega \mu_0)$
|
|
// After some derivation, we get:
|
|
// $F{S^-1} = \delta(t) - \sigma/\mu_0exp[-t(\sigma/\mu_0)]u(t)$
|
|
//```
|
|
let sigma = pseudo_conductivity;
|
|
|
|
let e = (-sigma * twice_delta_t).exp();
|
|
let aux = &mut pml.aux_h;
|
|
// (S^-1 \nabla) x E
|
|
Vec3::from((
|
|
aux.sy_conv_dfz_dy.step_only_sigma(e.y(), delta_ez_y)
|
|
- aux.sz_conv_dfy_dz.step_only_sigma(e.z(), delta_ey_z),
|
|
aux.sz_conv_dfx_dz.step_only_sigma(e.z(), delta_ex_z)
|
|
- aux.sx_conv_dfz_dx.step_only_sigma(e.x(), delta_ez_x),
|
|
aux.sx_conv_dfy_dx.step_only_sigma(e.x(), delta_ey_x)
|
|
- aux.sy_conv_dfx_dy.step_only_sigma(e.y(), delta_ex_y),
|
|
))
|
|
}
|
|
};
|
|
|
|
// From: $\nabla x E = -dB/dt$
|
|
let db_dt = -nabla_e;
|
|
// 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 state = CellState {
|
|
e: *self.e,
|
|
h: *self.h,
|
|
};
|
|
self.mat.step_b(&state, delta_b);
|
|
// mu0 (H + M) = B
|
|
// therefore, H = B mu0^-1 - M
|
|
// println!("cpu-step_h: nabla_e: {}", nabla_e);
|
|
// println!("cpu-step_h: delta_h: {}", delta_b * R::mu0_inv());
|
|
*self.h = (b_prev + delta_b) * R::mu0_inv() - self.mat.m();
|
|
}
|
|
}
|
|
|
|
impl<'a, R: Real, M: Material<R>> ECell<'a, R, M> {
|
|
/// delta_t = timestep covered by this function. i.e. it should be half the timestep of the simulation
|
|
/// feature_size = how many units apart is the center of each adjacent cell on the grid.
|
|
fn step(&mut self, neighbors: Neighbors<'_, R>, delta_t: R, inv_feature_size: R) {
|
|
// ```tex
|
|
// Ampere's circuital law with Maxwell's addition, in SI units ("macroscopic version"):
|
|
// $\nabla x H = J_f + dD/dt$ where $J_f$ = current density = $\sigma E$, $\sigma$ being a material parameter ("conductivity")
|
|
// Note that $D = \epsilon_0 E + P$, but we don't simulate any material where $P \ne 0$.
|
|
// Expand $D$: $\nabla x H = J_f + \epsilon_0 dE/dt$
|
|
// Expand $J$: $\nabla x H = \sigma E + \epsilon_0 dE/dt$
|
|
// Rearrange: $dE/dt = 1/\epsilon_0 (\nabla x H - \sigma E)$
|
|
//
|
|
// let $E_p$ be $E$ at $T-\Delta t$ and $E_n$ be $E$ at $T+\Delta t$.
|
|
// $(E_n - E_p)/(2\Delta{t}) = 1/\epsilon_0 (\nabla x H - \sigma (E_n + E_p)/2)$
|
|
// Normalize: $E_n - E_p = 2\Delta{t}/\epsilon_0 (\nabla x H - \sigma (E_n + E_p)/2)$
|
|
// Expand: $E_n - E_p = 2\Delta{t}/\epsilon_0 \nabla x H - \sigma \Delta{t}/\epsilon_0 (E_n + E_p)$
|
|
// Rearrange: $E_n(1 + \sigma \Delta{t}/\epsilon_0) = E_p(1 - \sigma \Delta{t}/\epsilon_0) + 2\Delta{t}/\epsilon_0 \nabla x H$
|
|
// Then $E_n$ (i.e. the E vector after this step) is trivially solved
|
|
//
|
|
// Note: $\nabla x H$ evaluates to:
|
|
// | i j k | [ dHz/dy - dHy/dz ]
|
|
// | | | |
|
|
// | d/dx d/dy d/dz | = | dHx/dz - dHz/dx |
|
|
// | | | |
|
|
// | Hx Hy Hz | [ dHy/dx - dHx/dy ]
|
|
//
|
|
// Boundary conditions (TODO):
|
|
// For something on the yz plane, E(x=0, t+1) should be E(x=1, t). i.e. the wave travels
|
|
// $\Delta x$ = 1 per timestep. Roughly, we want to take the line from E(x=0, t), E(x=1, t)
|
|
// and continue that to E(x=-1, t) ?
|
|
//
|
|
// ```
|
|
|
|
let twice_delta_t = R::two() * delta_t;
|
|
|
|
let e_prev = *self.e;
|
|
|
|
let (delta_hy_x, delta_hz_x) = match (neighbors.left, neighbors.right) {
|
|
(Some(left), _) => (
|
|
inv_feature_size * (self.h.y() - left.y()),
|
|
inv_feature_size * (self.h.z() - left.z()),
|
|
),
|
|
_ => (R::zero(), R::zero()),
|
|
};
|
|
let (delta_hx_y, delta_hz_y) = match (neighbors.up, neighbors.down) {
|
|
(Some(up), _) => (
|
|
inv_feature_size * (self.h.x() - up.x()),
|
|
inv_feature_size * (self.h.z() - up.z()),
|
|
),
|
|
_ => (R::zero(), R::zero()),
|
|
};
|
|
let (delta_hx_z, delta_hy_z) = match (neighbors.out, neighbors.in_) {
|
|
(Some(out), _) => (
|
|
inv_feature_size * (self.h.x() - out.x()),
|
|
inv_feature_size * (self.h.y() - out.y()),
|
|
),
|
|
_ => (R::zero(), R::zero()),
|
|
};
|
|
|
|
let StepParametersMut { conductivity, pml } = self.mat.step_parameters_mut();
|
|
let nabla_h = match pml {
|
|
None => {
|
|
// \nabla x H
|
|
Vec3::from((delta_hz_y - delta_hy_z, delta_hx_z - delta_hz_x, delta_hy_x - delta_hx_y))
|
|
}
|
|
Some((pml, PmlParameters { pseudo_conductivity })) => {
|
|
//```tex
|
|
// S = propagation_const + pseudo_conductivity / (dc_shift + j\omega \epsilon_0)
|
|
// $S = 1 + sigma/(j\omega \epsilon_0)$
|
|
// After some derivation, we get:
|
|
// $F{S^-1} = \delta(t) - \sigma/\epsilon_0exp[-t(\sigma/\epsilon_0)]u(t)$
|
|
//```
|
|
let sigma = pseudo_conductivity;
|
|
|
|
let e = (-sigma * twice_delta_t).exp();
|
|
let aux = &mut pml.aux_e;
|
|
// (S^-1 \nabla) x H
|
|
Vec3::from((
|
|
aux.sy_conv_dfz_dy.step_only_sigma(e.y(), delta_hz_y)
|
|
- aux.sz_conv_dfy_dz.step_only_sigma(e.z(), delta_hy_z),
|
|
aux.sz_conv_dfx_dz.step_only_sigma(e.z(), delta_hx_z)
|
|
- aux.sx_conv_dfz_dx.step_only_sigma(e.x(), delta_hz_x),
|
|
aux.sx_conv_dfy_dx.step_only_sigma(e.x(), delta_hy_x)
|
|
- aux.sy_conv_dfx_dy.step_only_sigma(e.y(), delta_hx_y),
|
|
))
|
|
}
|
|
};
|
|
|
|
// From: $\nabla x H = \epsilon_0 dE/dt + \sigma E$
|
|
let de_dt = solve_step_diff_eq_linear(
|
|
nabla_h,
|
|
Vec3::uniform(R::eps0()), // dE/dt coeff
|
|
conductivity, // E coeff
|
|
delta_t,
|
|
e_prev,
|
|
);
|
|
|
|
// println!("cpu-step_e: nabla_h: {}", nabla_h);
|
|
// println!("cpu-step_e: delta_e: {}", de_dt*twice_delta_t);
|
|
*self.e = e_prev + de_dt*twice_delta_t;
|
|
}
|
|
}
|
|
|
|
#[derive(Copy, Clone, Debug, Default, PartialEq, Serialize, Deserialize)]
|
|
pub struct CellState<R> {
|
|
e: Vec3<R>,
|
|
h: Vec3<R>,
|
|
}
|
|
|
|
impl<R: Real> CellState<R> {
|
|
pub fn e(&self) -> Vec3<R> {
|
|
self.e
|
|
}
|
|
pub fn h(&self) -> Vec3<R> {
|
|
self.h
|
|
}
|
|
|
|
pub fn with_m(&self, m: Vec3<R>) -> CellStateWithM<R> {
|
|
CellStateWithM {
|
|
e: self.e,
|
|
h: self.h,
|
|
m,
|
|
}
|
|
}
|
|
}
|
|
|
|
#[derive(Copy, Clone, Debug, Default, PartialEq, Serialize, Deserialize)]
|
|
pub struct CellStateWithM<R> {
|
|
e: Vec3<R>,
|
|
h: Vec3<R>,
|
|
m: Vec3<R>,
|
|
}
|
|
|
|
impl<R: Real> CellStateWithM<R> {
|
|
pub fn e(&self) -> Vec3<R> {
|
|
self.e
|
|
}
|
|
pub fn h(&self) -> Vec3<R> {
|
|
self.h
|
|
}
|
|
pub fn m(&self) -> Vec3<R> {
|
|
self.m
|
|
}
|
|
pub fn b(&self) -> Vec3<R> {
|
|
(self.h() + self.m()) * R::mu0()
|
|
}
|
|
}
|
|
|
|
/// Most of these serve as "smoke tests", to catch and pinpoint the most egregious of errors, or to
|
|
/// understand the impact of future changes (like changing default boundary conditions).
|
|
#[cfg(test)]
|
|
mod test {
|
|
use super::*;
|
|
use crate::geom::{Cube, Region, WorldRegion};
|
|
use crate::mat::{Conductor, Pml, Static};
|
|
use crate::real::{R64, ToFloat as _};
|
|
use float_eq::assert_float_eq;
|
|
use more_asserts::*;
|
|
use rand::rngs::StdRng;
|
|
use rand::{Rng as _, SeedableRng as _};
|
|
use std::fmt::Debug;
|
|
use std::ops::RangeBounds;
|
|
|
|
fn do_conv_flt_test(a: f64, b: f64, c: f64, delta_t: f64, signal: &[f64], expected: &[f64]) {
|
|
let mut cf = ConvFlt(0.0);
|
|
for (s, e) in signal.iter().zip(expected.iter()) {
|
|
let v = cf.step(
|
|
a,
|
|
b,
|
|
c,
|
|
*s,
|
|
delta_t,
|
|
);
|
|
assert_float_eq!(v, e, r2nd <= 1e-6);
|
|
}
|
|
}
|
|
|
|
#[test]
|
|
fn conv_flt_trivial() {
|
|
let signal = [2.0, 7.0, -3.0];
|
|
// kernel: 1.0 \delta(t)
|
|
let expected = [2.0, 7.0, -3.0];
|
|
let (a, b, c, delta_t) = (1.0, 0.0, 1.0, 0.1);
|
|
do_conv_flt_test(a, b, c, delta_t, &signal, &expected);
|
|
}
|
|
|
|
#[test]
|
|
fn conv_flt_trivial_scaled() {
|
|
let signal = [2.0, 7.0, -3.0];
|
|
// kernel: 10.0 \delta(t)
|
|
let expected = [20.0, 70.0, -30.0];
|
|
let (a, b, c, delta_t) = (10.0, 0.0, 1.0, 0.1);
|
|
do_conv_flt_test(a, b, c, delta_t, &signal, &expected);
|
|
}
|
|
|
|
#[test]
|
|
fn conv_flt_exp_trivial() {
|
|
let signal = [2.0, 0.0, 0.0];
|
|
// kernel: e(-t)
|
|
// \int_0^1 e(-t) dt = [1 - exp(-1)]
|
|
let exp_neg_0 = (-0.0f64).exp();
|
|
let exp_neg_1 = (-1.0f64).exp();
|
|
let exp_neg_2 = (-2.0f64).exp();
|
|
let exp_neg_3 = (-3.0f64).exp();
|
|
let expected = [
|
|
2.0*(exp_neg_0 - exp_neg_1),
|
|
2.0*(exp_neg_1 - exp_neg_2),
|
|
2.0*(exp_neg_2 - exp_neg_3),
|
|
];
|
|
let (a, b, c, delta_t) = (0.0, 1.0, 1.0, 1.0);
|
|
do_conv_flt_test(a, b, c, delta_t, &signal, &expected);
|
|
}
|
|
|
|
#[test]
|
|
fn conv_flt_exp_time_scaled() {
|
|
let signal = [2.0, 0.0, 0.0];
|
|
// kernel: e(-3*t)
|
|
// \int_0^0.2 e(-3*t) dt = [1 - exp(-0.6)]/3
|
|
let exp_neg_00 = (-0.0f64).exp();
|
|
let exp_neg_06 = (-0.6f64).exp();
|
|
let exp_neg_12 = (-1.2f64).exp();
|
|
let exp_neg_18 = (-1.8f64).exp();
|
|
let expected = [
|
|
2.0/3.0*(exp_neg_00 - exp_neg_06),
|
|
2.0/3.0*(exp_neg_06 - exp_neg_12),
|
|
2.0/3.0*(exp_neg_12 - exp_neg_18),
|
|
];
|
|
let (a, b, c, delta_t) = (0.0, 1.0, 3.0, 0.2);
|
|
do_conv_flt_test(a, b, c, delta_t, &signal, &expected);
|
|
}
|
|
|
|
#[test]
|
|
fn conv_flt_exp_multi_input() {
|
|
let signal = [2.0, 7.0, -3.0];
|
|
// kernel: e(-3*t)
|
|
// \int_0^0.2 e(-3*t) dt = [1 - exp(-0.6)]/3
|
|
let exp_neg_00 = (-0.0f64).exp();
|
|
let exp_neg_06 = (-0.6f64).exp();
|
|
let exp_neg_12 = (-1.2f64).exp();
|
|
let exp_neg_18 = (-1.8f64).exp();
|
|
let expected = [
|
|
2.0/3.0*(exp_neg_00 - exp_neg_06),
|
|
7.0/3.0*(exp_neg_00 - exp_neg_06) + 2.0/3.0*(exp_neg_06 - exp_neg_12),
|
|
-3.0/3.0*(exp_neg_00 - exp_neg_06)
|
|
+ 7.0/3.0*(exp_neg_06 - exp_neg_12)
|
|
+ 2.0/3.0*(exp_neg_12 - exp_neg_18),
|
|
];
|
|
let (a, b, c, delta_t) = (0.0, 1.0, 3.0, 0.2);
|
|
do_conv_flt_test(a, b, c, delta_t, &signal, &expected);
|
|
}
|
|
|
|
#[test]
|
|
fn conv_flt_all_params() {
|
|
let signal = [2.0, 7.0, -3.0];
|
|
// kernel: 0.3 \delta(t) + 0.4 * e(-3*t)
|
|
// \int_0^0.2 e(-3*t) dt = [1 - exp(-0.6)]/3
|
|
let exp_neg_00 = (-0.0f64).exp();
|
|
let exp_neg_06 = (-0.6f64).exp();
|
|
let exp_neg_12 = (-1.2f64).exp();
|
|
let exp_neg_18 = (-1.8f64).exp();
|
|
let expected_exp = [
|
|
2.0/3.0*(exp_neg_00 - exp_neg_06),
|
|
7.0/3.0*(exp_neg_00 - exp_neg_06) + 2.0/3.0*(exp_neg_06 - exp_neg_12),
|
|
-3.0/3.0*(exp_neg_00 - exp_neg_06)
|
|
+ 7.0/3.0*(exp_neg_06 - exp_neg_12)
|
|
+ 2.0/3.0*(exp_neg_12 - exp_neg_18),
|
|
];
|
|
let expected = [
|
|
0.3 * signal[0] + 0.4 * expected_exp[0],
|
|
0.3 * signal[1] + 0.4 * expected_exp[1],
|
|
0.3 * signal[2] + 0.4 * expected_exp[2],
|
|
];
|
|
let (a, b, c, delta_t) = (0.3, 0.4, 3.0, 0.2);
|
|
do_conv_flt_test(a, b, c, delta_t, &signal, &expected);
|
|
}
|
|
|
|
fn assert_vec3_eq<R: Real>(actual: Vec3<R>, expected: Vec3<R>) {
|
|
let actual = actual.cast::<f64>();
|
|
let expected = expected.cast::<f64>();
|
|
assert_float_eq!(actual.x(), expected.x(), r2nd <= 1e-6);
|
|
assert_float_eq!(actual.y(), expected.y(), r2nd <= 1e-6);
|
|
assert_float_eq!(actual.z(), expected.z(), r2nd <= 1e-6);
|
|
}
|
|
|
|
#[test]
|
|
fn step_diff_eq_trivial() {
|
|
// f(t) = 0.0
|
|
// nabla_g = 0.1 df/dt
|
|
let df_dt = solve_step_diff_eq(
|
|
Vec3::new(2.0, 3.0, 4.0), // nabla_g
|
|
Vec3::uniform(0.1), // df_dt_coeff
|
|
Vec3::zero(),
|
|
Vec3::zero(),
|
|
Vec3::zero(),
|
|
0.5, // delta_t
|
|
Vec3::zero(),
|
|
Vec3::zero(),
|
|
Vec3::zero(),
|
|
);
|
|
assert_vec3_eq(df_dt, Vec3::new(20.0, 30.0, 40.0));
|
|
}
|
|
|
|
#[test]
|
|
fn step_diff_eq_trivial_3rd_order() {
|
|
// f(t) = 0.0
|
|
// nabla_g = 0.1 df/dt - f - 3 f_int + 0.4 f_int_int
|
|
let df_dt = solve_step_diff_eq(
|
|
Vec3::new(2.0, 3.0, 4.0), // nabla_g
|
|
Vec3::uniform(0.1), // df_dt_coeff
|
|
Vec3::uniform(-1.0), // f_coeff
|
|
Vec3::uniform(-3.0), // f_int_coeff
|
|
Vec3::uniform(0.4), // f_int_int_coeff
|
|
0.0, // delta_t
|
|
Vec3::zero(),
|
|
Vec3::zero(),
|
|
Vec3::zero(),
|
|
);
|
|
assert_vec3_eq(df_dt, Vec3::new(20.0, 30.0, 40.0));
|
|
}
|
|
|
|
#[test]
|
|
fn step_diff_eq_1st_order_no_dt() {
|
|
// nabla_g = 0.1 df/dt + 2 f
|
|
// Let f(t) = sin(2 t)
|
|
// then f'(t) = 2 cos(2 t)
|
|
let f = Vec3::uniform((5.0f64).sin());
|
|
let df_dt = Vec3::uniform(2.0 * (5.0f64).cos());
|
|
let nabla_g = df_dt*0.1 + f*2.0;
|
|
let actual_df_dt = solve_step_diff_eq(
|
|
nabla_g,
|
|
Vec3::uniform(0.1),
|
|
Vec3::uniform(2.0),
|
|
Vec3::zero(),
|
|
Vec3::zero(),
|
|
0.0, // delta_t
|
|
f, // f_prev
|
|
Vec3::zero(),
|
|
Vec3::zero(),
|
|
);
|
|
assert_vec3_eq(actual_df_dt, df_dt);
|
|
}
|
|
|
|
#[test]
|
|
fn step_diff_eq_2nd_order_no_dt() {
|
|
// nabla_g = 0.1 df/dt + 2 f + 0.4 f_int
|
|
// Let f(t) = sin(2 t)
|
|
// then f'(t) = 2 cos(2 t)
|
|
// then \int f(t) = -1/2 cos(2 t)
|
|
let f = Vec3::uniform((5.0f64).sin());
|
|
let df_dt = Vec3::uniform(2.0 * (5.0f64).cos());
|
|
let f_int = Vec3::uniform(-0.5 * (5.0f64).cos());
|
|
let nabla_g = df_dt*0.1 + f*2.0 + f_int*0.4;
|
|
let actual_df_dt = solve_step_diff_eq(
|
|
nabla_g,
|
|
Vec3::uniform(0.1),
|
|
Vec3::uniform(2.0),
|
|
Vec3::uniform(0.4),
|
|
Vec3::zero(),
|
|
0.0, // delta_t
|
|
f, // f_prev
|
|
f_int, // f_int_prev
|
|
Vec3::zero(),
|
|
);
|
|
assert_vec3_eq(actual_df_dt, df_dt);
|
|
}
|
|
|
|
#[test]
|
|
fn step_diff_eq_3rd_order_no_dt() {
|
|
// nabla_g = 0.1 df/dt + 2 f + 0.4 f_int + 0.3 f_int_int
|
|
// Let f(t) = sin(2 t)
|
|
// then f'(t) = 2 cos(2 t)
|
|
// then \int f(t) = -1/2 cos(2 t)
|
|
// then \int \int f(t) = -1/4 sin(2 t)
|
|
let f = Vec3::unit()*(5.0f64).sin();
|
|
let df_dt = Vec3::unit()*2.0 * (5.0f64).cos();
|
|
let f_int = Vec3::unit()*-0.5 * (5.0f64).cos();
|
|
let f_int_int = Vec3::unit()*-0.25 * (5.0f64).sin();
|
|
let nabla_g = df_dt*0.1 + f*2.0 + f_int*0.4 + f_int_int*0.3;
|
|
let actual_df_dt = solve_step_diff_eq(
|
|
nabla_g,
|
|
Vec3::uniform(0.1),
|
|
Vec3::uniform(2.0),
|
|
Vec3::uniform(0.4),
|
|
Vec3::uniform(0.3),
|
|
0.0, // delta_t
|
|
f, // f_prev
|
|
f_int, // f_int_prev
|
|
f_int_int, // f_int_int_prev
|
|
);
|
|
assert_vec3_eq(actual_df_dt, df_dt);
|
|
}
|
|
|
|
#[test]
|
|
fn step_diff_eq_1st_order_homogenous() {
|
|
// 0 = f'(t) + f(t)
|
|
// solution: f(t) = k*e^-t
|
|
let f = |t: f32| -> Vec3<f32> {
|
|
Vec3::unit()*(-t).exp()
|
|
};
|
|
let df_dt = |t: f32| -> Vec3<f32> {
|
|
Vec3::unit()*-(-t).exp()
|
|
};
|
|
let f_int = |t: f32| -> Vec3<f32> {
|
|
Vec3::unit()*-(-t).exp()
|
|
};
|
|
let f_int_int = |t: f32| -> Vec3<f32> {
|
|
Vec3::unit()*(-t).exp()
|
|
};
|
|
let tl = 3.0000;
|
|
let tm = 3.0001;
|
|
let actual_df_dt = solve_step_diff_eq(
|
|
Vec3::zero(),
|
|
Vec3::unit(), // df_dt_coeff
|
|
Vec3::unit(), // f_coeff
|
|
Vec3::zero(),
|
|
Vec3::zero(),
|
|
tm - tl, // delta_t
|
|
f(tl), // f_prev
|
|
f_int(tl),
|
|
f_int_int(tl),
|
|
);
|
|
assert_vec3_eq(actual_df_dt, df_dt(tm));
|
|
}
|
|
|
|
#[test]
|
|
fn step_diff_eq_1st_order_nonhomogenous() {
|
|
// 1 = f'(t) + f(t)
|
|
// let f(t) = 1 + e^-t
|
|
let f = |t: f32| -> Vec3<f32> {
|
|
Vec3::unit() + Vec3::unit()*(-t).exp()
|
|
};
|
|
let df_dt = |t: f32| -> Vec3<f32> {
|
|
Vec3::unit()*-(-t).exp()
|
|
};
|
|
let f_int = |t: f32| -> Vec3<f32> {
|
|
Vec3::unit()*-(-t).exp()
|
|
};
|
|
let f_int_int = |t: f32| -> Vec3<f32> {
|
|
Vec3::unit()*(-t).exp()
|
|
};
|
|
let tl = 3.0000;
|
|
let tm = 3.0001;
|
|
let actual_df_dt = solve_step_diff_eq(
|
|
Vec3::unit(), // nabla_g
|
|
Vec3::unit(), // df_dt_coeff
|
|
Vec3::unit(), // f_coeff
|
|
Vec3::zero(),
|
|
Vec3::zero(),
|
|
tm - tl, // delta_t
|
|
f(tl), // f_prev
|
|
f_int(tl),
|
|
f_int_int(tl),
|
|
);
|
|
assert_vec3_eq(actual_df_dt, df_dt(tm));
|
|
}
|
|
|
|
#[test]
|
|
fn step_diff_eq_3rd_order_nonhomogenous() {
|
|
// c0 = af'(t) + bf(t) + c\int f(t) + d\int \int f(t)
|
|
// let f(t) = v0*e^wt, and \int \int f(t) absorbs C
|
|
// These constants are chosen arbitrarily: this test
|
|
// should pass for any "valid" values
|
|
let v0 = Vec3::new(3.0, 5.0, -0.5);
|
|
let w = Vec3::new(0.1, -2.0, 0.4);
|
|
let c0 = Vec3::new(-1.0, -2.0, -3.0);
|
|
let a = Vec3::new(0.1, 1.0, 2.0); // must be non-zero
|
|
let b = Vec3::new(0.5, -0.3, 0.0);
|
|
let c = Vec3::new(-5.0, 0.0, -0.1);
|
|
let d = -w.elem_mul(c + w.elem_mul(b + w.elem_mul(a)));
|
|
let f = |t: f64| -> Vec3<f64> {
|
|
v0.elem_mul((w*t).exp())
|
|
};
|
|
let df_dt = |t: f64| -> Vec3<f64> {
|
|
v0.elem_mul((w*t).exp()).elem_mul(w)
|
|
};
|
|
let f_int = |t: f64| -> Vec3<f64> {
|
|
v0.elem_mul((w*t).exp()).elem_div(w)
|
|
};
|
|
let f_int_int = |t: f64| -> Vec3<f64> {
|
|
v0.elem_mul((w*t).exp()).elem_div(w).elem_div(w) + c0.elem_div(d)
|
|
};
|
|
let tl = 3.00000000;
|
|
let tm = 3.00000001;
|
|
let actual_df_dt = solve_step_diff_eq(
|
|
c0, // nabla_g
|
|
a, // df_dt_coeff
|
|
b, // f_coeff
|
|
c, // int_f_coeff
|
|
d, // int_int_f_coeff
|
|
tm - tl, // delta_t
|
|
f(tl), // f_prev
|
|
f_int(tl), // f_int_prev
|
|
f_int_int(tl), // f_int_int_prev
|
|
);
|
|
assert_vec3_eq(actual_df_dt, df_dt(tm));
|
|
}
|
|
|
|
fn energy(s: &dyn SampleableSim) -> f32 {
|
|
0.5 * s.feature_volume() * s.map_sum_enumerated(|_pos: Index, cell| {
|
|
// println!("{}: {}", _pos, cell.e());
|
|
cell.e().mag_sq().to_f64()
|
|
}).to_f32()
|
|
}
|
|
|
|
fn energy_now_and_then<S: GenericSim>(state: &mut S, frames: u32) -> (f32, f32) {
|
|
let energy_0 = energy(state);
|
|
for _f in 0..frames {
|
|
// println!("e({}) = {}", _f, energy(state));
|
|
state.step();
|
|
}
|
|
let energy_1 = energy(state);
|
|
(energy_0, energy_1)
|
|
}
|
|
|
|
#[test]
|
|
fn accessors() {
|
|
let mut state = SimState::<R64, Static<R64>>::new(Index((15, 17, 19).into()), 1e-6);
|
|
assert_eq!(state.width(), 15);
|
|
assert_eq!(state.height(), 17);
|
|
assert_eq!(state.depth(), 19);
|
|
assert_eq!(state.time(), 0.0);
|
|
// feat_size / C * courant: 1e-6 / 3e8 * 0.577
|
|
assert_eq!(state.timestep(), 0.000000000000001924664829293337);
|
|
state.step();
|
|
assert_eq!(state.time(), state.timestep());
|
|
}
|
|
|
|
#[test]
|
|
fn energy_conservation_over_time() {
|
|
let mut state = SimState::<R64, Static<R64>>::new(Index((2001, 1, 1).into()), 1e-6);
|
|
for t in 0..100 {
|
|
let angle = (t as f32)/100.0*f32::two_pi();
|
|
let sig = angle.sin();
|
|
let gate = 0.5*(1.0 - angle.cos());
|
|
//println!("{}", g.exp());
|
|
state.impulse_ez(Index((1000, 0, 0).into()), gate*sig);
|
|
//dyn_state.impulse_by(index((x, 0, 0).into()), 1.5e-9*g.exp());
|
|
state.step();
|
|
}
|
|
let (energy_0, energy_1) = energy_now_and_then(&mut state, 800);
|
|
// This has some sensitivity to courant number. But surprisingly, not much to f32 v.s. f64
|
|
assert_float_eq!(energy_1, energy_0, r2nd <= 1e-9);
|
|
}
|
|
|
|
#[test]
|
|
fn sane_boundary_conditions() {
|
|
let mut state = SimState::<R64, Static<R64>>::new(Index((21, 21, 21).into()), 1e-6);
|
|
for t in 0..40 {
|
|
let angle = (t as f32)/10.0*f32::two_pi();
|
|
let sig = angle.sin();
|
|
let gate = 1e9 * 0.5*(1.0 - angle.cos());
|
|
state.impulse_ez(Index((10, 10, 10).into()), gate*sig);
|
|
state.step();
|
|
}
|
|
let (energy_0, energy_1) = energy_now_and_then(&mut state, 500);
|
|
// Default boundary conditions reflect waves.
|
|
assert_float_eq!(energy_1, energy_0, r2nd <= 1e-2);
|
|
}
|
|
|
|
/// Fill the world with the provided material and a stimulus.
|
|
/// Measure energy at the start, and then again after advancing many steps.
|
|
/// Return these two measurements (energy(t=0), energy(t=~=1000))
|
|
fn conductor_test<M: Into<mat::Static<R64>>>(mat: M) -> (f32, f32) {
|
|
let mut state = SimState::<R64, Static<R64>>::new(Index((201, 1, 1).into()), 1e-6);
|
|
state.fill_region(&WorldRegion, mat.into());
|
|
for t in 0..100 {
|
|
let progress = (t as f32)/100.0;
|
|
let angle = 2.0*progress*f32::two_pi();
|
|
let sig = angle.sin();
|
|
let gate = (progress*f32::pi()).sin();
|
|
state.impulse_ez(Index((100, 0, 0).into()), gate*sig);
|
|
state.step();
|
|
}
|
|
energy_now_and_then(&mut state, 1000)
|
|
}
|
|
|
|
#[test]
|
|
fn conductor_dissipates_energy() {
|
|
let (energy_0, energy_1) = conductor_test(Conductor::new(R64::from_f32(1e3)));
|
|
assert_float_eq!(energy_1/energy_0, 0.0, abs <= 1e-6);
|
|
}
|
|
|
|
#[test]
|
|
fn anisotropic_conductor_inactive_x() {
|
|
let (energy_0, energy_1) = conductor_test(Conductor::new_anisotropic(
|
|
Vec3::new_x(1e3)
|
|
));
|
|
assert_gt!(energy_1, 0.9*energy_0);
|
|
// XXX: I think this gains energy because we only set the E field and not the H field
|
|
assert_lt!(energy_1, 1.5*energy_0);
|
|
}
|
|
|
|
#[test]
|
|
fn anisotropic_conductor_inactive_y() {
|
|
let (energy_0, energy_1) = conductor_test(Conductor::new_anisotropic(
|
|
Vec3::new_y(1e3)
|
|
));
|
|
assert_gt!(energy_1, 0.9*energy_0);
|
|
assert_lt!(energy_1, 1.5*energy_0);
|
|
}
|
|
|
|
#[test]
|
|
fn anisotropic_conductor_active_z() {
|
|
let (energy_0, energy_1) = conductor_test(Conductor::new_anisotropic(
|
|
Vec3::new_z(1e3)
|
|
));
|
|
assert_float_eq!(energy_1/energy_0, 0.0, abs <= 1e-6);
|
|
}
|
|
|
|
fn state_for_pml(size: Index) -> SimState<R64, Pml<R64>> {
|
|
let mut state = SimState::new(size, 1e-6);
|
|
let timestep = state.timestep();
|
|
state.fill_boundary_using(size/4, |boundary_ness| {
|
|
let b = boundary_ness.elem_pow(3.0);
|
|
let conductivity = b * (0.5 / timestep);
|
|
// K scales up to 11; dc_shift scales from 0.05:
|
|
// https://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.461.4606&rep=rep1&type=pdf
|
|
// let propagation = Vec3::unit() + b * 10.0;
|
|
// let dc_shift = (Vec3::unit() - b) * 0.05;
|
|
//
|
|
// let propagation = Vec3::unit();
|
|
// let dc_shift = Vec3::unit() * 1e-6;
|
|
Pml::new(conductivity)
|
|
});
|
|
state
|
|
}
|
|
|
|
/// Like state_for_pml, but the boundary is an ordinary electric conductor -- not PML
|
|
fn state_for_graded_conductor(size: Index) -> SimState<R64, Conductor<R64>> {
|
|
let mut state = SimState::new(size, 1e-6);
|
|
let timestep = state.timestep();
|
|
state.fill_boundary_using(size/4, |boundary_ness| {
|
|
let b = boundary_ness.elem_pow(3.0);
|
|
let conductivity = f32::eps0() * b.mag() * 0.5 / timestep;
|
|
Conductor::new(conductivity.cast())
|
|
});
|
|
state
|
|
}
|
|
|
|
/// Stimulus which applies a (possibly different) e*sin(wt)*gate(t) at every point
|
|
struct PmlStim<F> {
|
|
/// Maps location -> (field direction, frequency)
|
|
f: F,
|
|
t_end: f32,
|
|
feat_size: f32,
|
|
sqrt_vol: f32,
|
|
}
|
|
impl<F: Fn(Index) -> (Vec3<f32>, f32) + Sync> AbstractStimulus for PmlStim<F> {
|
|
fn at(&self, t_sec: f32, pos: Meters) -> (Vec3<f32>, Vec3<f32>) {
|
|
let angle = t_sec/self.t_end * f32::two_pi();
|
|
let gate = 0.5*(1.0 - angle.cos());
|
|
let (e, hz) = (self.f)(pos.to_index(self.feat_size));
|
|
let sig_angle = angle*hz;
|
|
let sig = sig_angle.sin();
|
|
// TODO: test H component?
|
|
(e*(gate*sig / (self.t_end * self.sqrt_vol)), Vec3::zero())
|
|
}
|
|
}
|
|
|
|
/// Returns the energy ratio (Eend / Estart)
|
|
fn pml_test_full_interior<F, S: GenericSim>(state: &mut S, f: F) -> f32
|
|
where F: Fn(Index) -> (Vec3<f32>, f32) + Sync // returns (E vector, cycles (like Hz))
|
|
{
|
|
let stim = PmlStim {
|
|
f,
|
|
t_end: 50.0 * state.timestep(),
|
|
feat_size: state.feature_size(),
|
|
sqrt_vol: state.volume().sqrt(),
|
|
};
|
|
for _t in 0..50 {
|
|
// println!("estim({}) = {}", _t, energy(state));
|
|
state.step_multiple(1, &stim);
|
|
}
|
|
let (energy_0, energy_1) = energy_now_and_then(state, 1000);
|
|
// sanity check the starting energy
|
|
assert_gt!(energy_0, 1e-11);
|
|
assert_lt!(energy_0, 1.0);
|
|
energy_1 / energy_0
|
|
}
|
|
|
|
fn pml_test_over_region<R, F, S>(state: &mut S, reg: R, f: F) -> f32
|
|
where
|
|
R: Region,
|
|
F: Fn(Index) -> (Vec3<f32>, f32) + Sync,
|
|
S: GenericSim
|
|
{
|
|
let feat = state.feature_size();
|
|
pml_test_full_interior(state, |idx| {
|
|
if reg.contains(idx.to_meters(feat)) {
|
|
f(idx)
|
|
} else {
|
|
(Vec3::zero(), 0.0)
|
|
}
|
|
})
|
|
}
|
|
fn pml_test_uniform_over_region<R: Region, S: GenericSim>(state: &mut S, reg: R, e: Vec3<f32>, hz: f32) -> f32 {
|
|
pml_test_over_region(state, reg, |_idx| {
|
|
(e, hz)
|
|
})
|
|
}
|
|
|
|
fn pml_test_at<S: GenericSim>(state: &mut S, e: Vec3<f32>, center: Index) -> f32 {
|
|
pml_test_full_interior(state, |idx| {
|
|
if idx == center {
|
|
(e, 1.0)
|
|
} else {
|
|
(Vec3::zero(), 1.0)
|
|
}
|
|
})
|
|
}
|
|
|
|
fn pml_test<S: GenericSim, R: RangeBounds<f32> + Debug>(state: &mut S, e: Vec3<f32>, en_range: R) {
|
|
let center = Index(state.size().0 / 2);
|
|
let pml = pml_test_at(state, e, center);
|
|
// PML should absorb all energy
|
|
assert!(en_range.contains(&pml), "{} not in {:?}", pml, en_range);
|
|
}
|
|
|
|
fn pml_test_against_baseline<P: GenericSim, B: GenericSim>(pml_state: &mut P, baseline_state: &mut B, e: Vec3<f32>) {
|
|
assert_eq!(pml_state.size(), baseline_state.size());
|
|
let center = Index(pml_state.size().0 / 2);
|
|
let en_pml = pml_test_at(pml_state, e, center);
|
|
let en_baseline = pml_test_at(baseline_state, e, center);
|
|
// PML should absorb all energy
|
|
println!("en_pml: {}", en_pml);
|
|
println!("en_baseline: {}", en_baseline);
|
|
assert_float_eq!(en_pml, 0.0, abs <= 1.2e-6);
|
|
// XXX: why is PML energy sometimes MORE here?
|
|
assert_lt!(en_pml / en_baseline, 3.0, "{}/{}", en_pml, en_baseline);
|
|
// panic!("debugging");
|
|
}
|
|
|
|
#[test]
|
|
fn pml_boundary_conditions_x() {
|
|
let mut pml_state = state_for_pml(Index::new(201, 1, 1));
|
|
let mut baseline_state = state_for_graded_conductor(Index::new(201, 1, 1));
|
|
pml_test_against_baseline(&mut pml_state, &mut baseline_state, Vec3::unit_z());
|
|
}
|
|
|
|
#[test]
|
|
fn pml_boundary_conditions_y() {
|
|
let mut pml_state = state_for_pml(Index::new(1, 201, 1));
|
|
let mut baseline_state = state_for_graded_conductor(Index::new(1, 201, 1));
|
|
pml_test_against_baseline(&mut pml_state, &mut baseline_state, Vec3::unit_x());
|
|
}
|
|
|
|
#[test]
|
|
fn pml_boundary_conditions_z() {
|
|
let mut pml_state = state_for_pml(Index::new(1, 1, 201));
|
|
let mut baseline_state = state_for_graded_conductor(Index::new(1, 1, 201));
|
|
pml_test_against_baseline(&mut pml_state, &mut baseline_state, Vec3::unit_y());
|
|
}
|
|
|
|
#[test]
|
|
fn pml_boundary_conditions_plane_xy() {
|
|
let mut pml_state = state_for_pml(Index::new(201, 201, 1));
|
|
let mut baseline_state = state_for_graded_conductor(Index::new(201, 201, 1));
|
|
pml_test_against_baseline(&mut pml_state, &mut baseline_state, Vec3::unit_z());
|
|
}
|
|
|
|
#[test]
|
|
fn pml_boundary_conditions_plane_xz() {
|
|
let mut pml_state = state_for_pml(Index::new(201, 1, 201));
|
|
let mut baseline_state = state_for_graded_conductor(Index::new(201, 1, 201));
|
|
pml_test_against_baseline(&mut pml_state, &mut baseline_state, Vec3::unit_y());
|
|
}
|
|
|
|
#[test]
|
|
fn pml_boundary_conditions_plane_yz() {
|
|
let mut pml_state = state_for_pml(Index::new(1, 201, 201));
|
|
let mut baseline_state = state_for_graded_conductor(Index::new(1, 201, 201));
|
|
pml_test_against_baseline(&mut pml_state, &mut baseline_state, Vec3::unit_x());
|
|
}
|
|
|
|
fn state_for_monodirectional_pml(size: Index) -> SimState<R64, Pml<R64>> {
|
|
let mut state = SimState::new(size, 1e-6);
|
|
let timestep = state.timestep();
|
|
state.fill_boundary_using(size/4, |boundary_ness| {
|
|
let b = boundary_ness.elem_pow(3.0);
|
|
let coord_stretch = b * (0.5 / timestep);
|
|
// let coord_stretch = b * 0.5;
|
|
// Force the stretch to be only along one axis
|
|
let coord_stretch: Vec3<_> = match coord_stretch.to_tuple() {
|
|
(x, y, z) if x >= y && x >= z => (x, 0.0, 0.0),
|
|
(x, y, z) if y >= x && y >= z => (0.0, y, 0.0),
|
|
(x, y, z) if z >= x && z >= y => (0.0, 0.0, z),
|
|
_ => unreachable!(),
|
|
}.into();
|
|
Pml::new(coord_stretch)
|
|
});
|
|
state
|
|
}
|
|
|
|
/// Despite PML existing, verify that it doesn't interfere with a specific wave.
|
|
fn pml_ineffective_mono_linear_test<F: Fn(Vec3<f32>) -> Vec3<f32>>(
|
|
size: Index, e: Vec3<f32>, shuffle: F
|
|
) {
|
|
let mut state = SimState::<R64, Pml<R64>>::new(size, 1e-6);
|
|
let timestep = state.timestep();
|
|
state.fill_boundary_using(size/4, |boundary_ness| {
|
|
let b = boundary_ness.elem_pow(3.0);
|
|
let cs = b * (0.5 / timestep);
|
|
// let cs = b * 0.5;
|
|
// permute the stretching so that it shouldn't interfere with the wave
|
|
let coord_stretch = shuffle(cs);
|
|
Pml::<R64>::new(coord_stretch)
|
|
});
|
|
let center = Index(state.size().0 / 2);
|
|
let en = pml_test_at(&mut state, e, center);
|
|
// XXX: would be better to compare it to an empty sim?
|
|
assert_float_eq!(en, 1.0, abs <= 0.1);
|
|
}
|
|
|
|
#[test]
|
|
fn pml_boundary_3d_inneffective_monodirectional_linear() {
|
|
for (size, e) in vec![
|
|
(Index::new(201, 1, 1), Vec3::unit_y()),
|
|
(Index::new(201, 1, 1), Vec3::unit_z()),
|
|
(Index::new(1, 201, 1), Vec3::unit_x()),
|
|
(Index::new(1, 201, 1), Vec3::unit_z()),
|
|
(Index::new(1, 1, 201), Vec3::unit_x()),
|
|
(Index::new(1, 1, 201), Vec3::unit_y()),
|
|
] {
|
|
pml_ineffective_mono_linear_test(
|
|
size, e, |v| Vec3::new(v.y(), v.z(), v.x())
|
|
);
|
|
pml_ineffective_mono_linear_test(
|
|
size, e, |v| Vec3::new(v.z(), v.x(), v.y())
|
|
);
|
|
}
|
|
}
|
|
|
|
#[test]
|
|
fn pml_boundary_3d_monodirectional_plane_xy() {
|
|
//let mut state = state_for_monodirectional_pml(Index::unit()*41);
|
|
let mut state = state_for_monodirectional_pml(Index::new(201, 201, 1));
|
|
pml_test(&mut state, Vec3::unit_z(), ..7e-5);
|
|
}
|
|
|
|
#[test]
|
|
fn pml_boundary_3d_monodirectional_plane_xz() {
|
|
//let mut state = state_for_monodirectional_pml(Index::unit()*41);
|
|
let mut state = state_for_monodirectional_pml(Index::new(201, 1, 201));
|
|
pml_test(&mut state, Vec3::unit_y(), ..7e-5);
|
|
}
|
|
|
|
#[test]
|
|
fn pml_boundary_3d_monodirectional_plane_yz() {
|
|
//let mut state = state_for_monodirectional_pml(Index::unit()*41);
|
|
let mut state = state_for_monodirectional_pml(Index::new(1, 201, 201));
|
|
pml_test(&mut state, Vec3::unit_x(), ..7e-5);
|
|
}
|
|
|
|
#[test]
|
|
fn pml_boundary_3d_monodirectional_x() {
|
|
let mut state = state_for_monodirectional_pml(Index::unit()*41);
|
|
pml_test(&mut state, Vec3::unit_x(), 1e30..);
|
|
}
|
|
|
|
#[test]
|
|
fn pml_boundary_3d_monodirectional_y() {
|
|
let mut state = state_for_monodirectional_pml(Index::unit()*41);
|
|
pml_test(&mut state, Vec3::unit_y(), 1e30..);
|
|
}
|
|
|
|
#[test]
|
|
fn pml_boundary_3d_monodirectional_z() {
|
|
let mut state = state_for_monodirectional_pml(Index::unit()*41);
|
|
pml_test(&mut state, Vec3::unit_z(), 1e30..);
|
|
}
|
|
|
|
#[test]
|
|
fn pml_boundary_3d_multidirectional() {
|
|
for dir in vec![Vec3::unit_x(), Vec3::unit_y(), Vec3::unit_z()] {
|
|
let mut pml_state = state_for_pml(Index::unit()*41);
|
|
let mut baseline_state = state_for_graded_conductor(Index::unit()*41);
|
|
pml_test_against_baseline(&mut pml_state, &mut baseline_state, dir);
|
|
}
|
|
}
|
|
|
|
#[test]
|
|
fn pml_boundary_3d_multidirectional_non_axis_aligned() {
|
|
for dir in vec![
|
|
Vec3::new(0.0, 1.0, 1.0),
|
|
Vec3::new(-1.0, 0.0, 1.0),
|
|
Vec3::new(1.0, -1.0, 0.0),
|
|
Vec3::new(1.0, -1.0, 1.0),
|
|
Vec3::new(-1.0, -1.0, -1.0),
|
|
Vec3::new(0.2, 0.5, 0.6),
|
|
Vec3::new(-2.0, 0.5, -0.5),
|
|
] {
|
|
let mut pml_state = state_for_pml(Index::unit()*41);
|
|
let mut baseline_state = state_for_graded_conductor(Index::unit()*41);
|
|
pml_test_against_baseline(&mut pml_state, &mut baseline_state, dir);
|
|
}
|
|
}
|
|
|
|
#[test]
|
|
fn pml_boundary_3d_multidirectional_off_center() {
|
|
for dir in vec![Vec3::unit_x(), Vec3::unit_y(), Vec3::unit_z()] {
|
|
for center in vec![
|
|
Index::new(15, 15, 15),
|
|
Index::new(15, 15, 25),
|
|
Index::new(15, 25, 15),
|
|
Index::new(15, 25, 25),
|
|
Index::new(25, 15, 15),
|
|
Index::new(25, 15, 25),
|
|
Index::new(25, 25, 15),
|
|
Index::new(25, 25, 25),
|
|
] {
|
|
let mut pml_state = state_for_pml(Index::unit()*41);
|
|
let mut baseline_state = state_for_graded_conductor(Index::unit()*41);
|
|
let pml_en = pml_test_at(&mut pml_state, dir, center);
|
|
let baseline_en = pml_test_at(&mut baseline_state, dir, center);
|
|
assert_float_eq!(pml_en, 0.0, abs <= 2e-5);
|
|
// XXX: why is PML energy MORE here?
|
|
assert_lt!(pml_en / baseline_en, 2.0);
|
|
}
|
|
}
|
|
}
|
|
|
|
// XXX the energies here are even worse than above.
|
|
// just don't use our PML implementation: it's not great.
|
|
// #[test]
|
|
// fn pml_boundary_3d_multidirectional_off_center_multisource() {
|
|
// for dir in vec![Vec3::unit_x(), Vec3::unit_y(), Vec3::unit_z()] {
|
|
// for cycles in vec![1.0, 2.0, 4.0, 8.0] {
|
|
// let mut pml_state = state_for_pml(Index::unit()*41);
|
|
// let mut baseline_state = state_for_graded_conductor(Index::unit()*41);
|
|
// let feat = pml_state.feature_size();
|
|
// let interior = Cube::new(Index::new(15, 15, 15).to_meters(feat), Index::new(25, 25, 25).to_meters(feat));
|
|
// let pml_en = pml_test_uniform_over_region(&mut pml_state, interior, dir, cycles);
|
|
// let baseline_en = pml_test_uniform_over_region(&mut baseline_state, interior, dir, cycles);
|
|
// assert_float_eq!(pml_en, 0.0, abs <= 4e-5);
|
|
// // XXX: why is PML energy MORE here?
|
|
// assert_lt!(pml_en / baseline_en, 5.0, "{} {} {}/{}", dir, cycles, pml_en, baseline_en);
|
|
// }
|
|
// }
|
|
// }
|
|
|
|
#[test]
|
|
fn pml_boundary_3d_chaotic_field() {
|
|
let excitation = |idx: Index| {
|
|
let seed = (idx.x() as u64) ^ ((idx.y() as u64) << 16) ^ ((idx.z() as u64) << 32);
|
|
let mut rng = StdRng::seed_from_u64(seed);
|
|
let dir = Vec3::new(
|
|
rng.gen_range(-1.0..1.0),
|
|
rng.gen_range(-1.0..1.0),
|
|
rng.gen_range(-1.0..1.0),
|
|
);
|
|
// XXX only works if it's a whole number. I suppose this makes sense though, as
|
|
// other numbers would create higher harmonics when gated.
|
|
let hz = rng.gen_range(0..=5);
|
|
(dir, hz as _)
|
|
};
|
|
let mut pml_state = state_for_pml(Index::unit()*41);
|
|
let mut baseline_state = state_for_graded_conductor(Index::unit()*41);
|
|
let feat = pml_state.feature_size();
|
|
let interior = Cube::new(Index::new(15, 15, 15).to_meters(feat), Index::new(25, 25, 25).to_meters(feat));
|
|
let en_pml = pml_test_over_region(&mut pml_state, interior, excitation);
|
|
let en_baseline = pml_test_over_region(&mut baseline_state, interior, excitation);
|
|
assert_float_eq!(en_pml, 0.0, abs <= 2e-5);
|
|
// XXX: why is PML energy MORE here?
|
|
assert_lt!(en_pml/en_baseline, 3.0);
|
|
}
|
|
}
|