rename 'coremem_types' -> 'coremem_cross' to better reflect its purpose
This commit is contained in:
382
crates/cross/src/step.rs
Normal file
382
crates/cross/src/step.rs
Normal file
@@ -0,0 +1,382 @@
|
||||
use core::ops::{Index, IndexMut};
|
||||
|
||||
use crate::compound::Optional;
|
||||
use crate::dim::DimensionedSlice;
|
||||
use crate::mat::Material;
|
||||
use crate::real::Real;
|
||||
use crate::vec::{Vec3, Vec3u};
|
||||
|
||||
#[cfg(feature = "serde")]
|
||||
use serde::{Serialize, Deserialize};
|
||||
|
||||
#[cfg_attr(feature = "serde", derive(Deserialize, Serialize))]
|
||||
#[cfg_attr(feature = "fmt", derive(Debug))]
|
||||
#[derive(Copy, Clone, Default)]
|
||||
pub struct SimMeta<R> {
|
||||
// TODO: make these private?
|
||||
pub dim: Vec3u,
|
||||
pub inv_feature_size: R,
|
||||
pub time_step: R,
|
||||
pub feature_size: R,
|
||||
}
|
||||
|
||||
impl<R: Copy> SimMeta<R> {
|
||||
pub fn dim(&self) -> Vec3u {
|
||||
self.dim
|
||||
}
|
||||
pub fn inv_feature_size(&self) -> R {
|
||||
self.inv_feature_size
|
||||
}
|
||||
pub fn time_step(&self) -> R {
|
||||
self.time_step
|
||||
}
|
||||
pub fn feature_size(&self) -> R {
|
||||
self.feature_size
|
||||
}
|
||||
}
|
||||
|
||||
impl<R: Real> SimMeta<R> {
|
||||
pub fn cast<R2: Real>(self) -> SimMeta<R2> {
|
||||
SimMeta {
|
||||
dim: self.dim,
|
||||
inv_feature_size: self.inv_feature_size.cast(),
|
||||
time_step: self.time_step.cast(),
|
||||
feature_size: self.feature_size.cast(),
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
/// Package the field vectors adjacent to some particular location.
|
||||
/// Particular those at negative offsets from the midpoint.
|
||||
/// This is used in step_e when looking at the H field deltas.
|
||||
#[derive(Copy, Clone)]
|
||||
pub struct VolumeSampleNeg<R> {
|
||||
pub mid: Vec3<R>,
|
||||
pub xm1: Optional<Vec3<R>>,
|
||||
pub ym1: Optional<Vec3<R>>,
|
||||
pub zm1: Optional<Vec3<R>>,
|
||||
}
|
||||
|
||||
impl<R: Copy + Default> VolumeSampleNeg<R> {
|
||||
pub fn from_indexable<I: Index<Vec3u, Output=Vec3<R>>>(i: &I, idx: Vec3u) -> Self {
|
||||
VolumeSampleNeg {
|
||||
mid: i[idx],
|
||||
xm1: prev_x(idx).map(|idx| i[idx]),
|
||||
ym1: prev_y(idx).map(|idx| i[idx]),
|
||||
zm1: prev_z(idx).map(|idx| i[idx]),
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
fn prev_x(idx: Vec3u) -> Optional<Vec3u> {
|
||||
match idx.into() {
|
||||
(0, _, _) => Optional::none(),
|
||||
(x, y, z) => Optional::some(Vec3u::new(x-1, y, z)),
|
||||
}
|
||||
}
|
||||
fn prev_y(idx: Vec3u) -> Optional<Vec3u> {
|
||||
match idx.into() {
|
||||
(_, 0, _) => Optional::none(),
|
||||
(x, y, z) => Optional::some(Vec3u::new(x, y-1, z)),
|
||||
}
|
||||
}
|
||||
fn prev_z(idx: Vec3u) -> Optional<Vec3u> {
|
||||
match idx.into() {
|
||||
(_, _, 0) => Optional::none(),
|
||||
(x, y, z) => Optional::some(Vec3u::new(x, y, z-1)),
|
||||
}
|
||||
}
|
||||
|
||||
impl<R: Real> VolumeSampleNeg<R> {
|
||||
/// Calculate the delta in H values amongst this cell and its neighbors (left/up/out)
|
||||
fn delta_h(self) -> FieldDeltas<R> {
|
||||
let mid = self.mid;
|
||||
// let (dfy_dx, dfz_dx) = self.xm1.map(|xm1| {
|
||||
// (mid.y() - xm1.y(), mid.z() - xm1.z())
|
||||
// }).unwrap_or_default();
|
||||
|
||||
// let (dfx_dy, dfz_dy) = self.ym1.map(|ym1| {
|
||||
// (mid.x() - ym1.x(), mid.z() - ym1.z())
|
||||
// }).unwrap_or_default();
|
||||
|
||||
// let (dfx_dz, dfy_dz) = self.zm1.map(|zm1| {
|
||||
// (mid.x() - zm1.x(), mid.y() - zm1.y())
|
||||
// }).unwrap_or_default();
|
||||
|
||||
let (dfy_dx, dfz_dx) = if self.xm1.is_some() {
|
||||
(mid.y() - self.xm1.unwrap().y(), mid.z() - self.xm1.unwrap().z())
|
||||
} else {
|
||||
(R::zero(), R::zero())
|
||||
};
|
||||
|
||||
let (dfx_dy, dfz_dy) = if self.ym1.is_some() {
|
||||
(mid.x() - self.ym1.unwrap().x(), mid.z() - self.ym1.unwrap().z())
|
||||
} else {
|
||||
(R::zero(), R::zero())
|
||||
};
|
||||
|
||||
let (dfx_dz, dfy_dz) = if self.zm1.is_some() {
|
||||
(mid.x() - self.zm1.unwrap().x(), mid.y() - self.zm1.unwrap().y())
|
||||
} else {
|
||||
(R::zero(), R::zero())
|
||||
};
|
||||
|
||||
FieldDeltas {
|
||||
dfy_dx,
|
||||
dfz_dx,
|
||||
dfx_dy,
|
||||
dfz_dy,
|
||||
dfx_dz,
|
||||
dfy_dz,
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
/// Package the field vectors adjacent to some particular location.
|
||||
/// Particular those at positive offsets from the midpoint.
|
||||
/// This is used in step_h when looking at the E field deltas.
|
||||
#[derive(Copy, Clone)]
|
||||
pub struct VolumeSamplePos<R> {
|
||||
pub mid: Vec3<R>,
|
||||
pub xp1: Optional<Vec3<R>>,
|
||||
pub yp1: Optional<Vec3<R>>,
|
||||
pub zp1: Optional<Vec3<R>>
|
||||
}
|
||||
|
||||
impl<R: Copy + Default> VolumeSamplePos<R> {
|
||||
pub fn from_indexable<I: Index<Vec3u, Output=Vec3<R>>>(i: &I, dim: Vec3u, idx: Vec3u) -> Self {
|
||||
VolumeSamplePos {
|
||||
mid: i[idx],
|
||||
xp1: next_x(dim, idx).map(|idx| i[idx]),
|
||||
yp1: next_y(dim, idx).map(|idx| i[idx]),
|
||||
zp1: next_z(dim, idx).map(|idx| i[idx]),
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
fn next_x(dim: Vec3u, idx: Vec3u) -> Optional<Vec3u> {
|
||||
match idx.into() {
|
||||
(x, y, z) if x + 1 < dim.x() => Optional::some(Vec3u::new(x+1, y, z)),
|
||||
_ => Optional::none(),
|
||||
}
|
||||
}
|
||||
fn next_y(dim: Vec3u, idx: Vec3u) -> Optional<Vec3u> {
|
||||
match idx.into() {
|
||||
(x, y, z) if y + 1 < dim.y() => Optional::some(Vec3u::new(x, y+1, z)),
|
||||
_ => Optional::none(),
|
||||
}
|
||||
}
|
||||
fn next_z(dim: Vec3u, idx: Vec3u) -> Optional<Vec3u> {
|
||||
match idx.into() {
|
||||
(x, y, z) if z + 1 < dim.z() => Optional::some(Vec3u::new(x, y, z+1)),
|
||||
_ => Optional::none(),
|
||||
}
|
||||
}
|
||||
|
||||
impl<R: Real> VolumeSamplePos<R> {
|
||||
/// Calculate the delta in E values amongst this cell and its neighbors (right/down/in)
|
||||
fn delta_e(self) -> FieldDeltas<R> {
|
||||
let mid = self.mid;
|
||||
// let (dfy_dx, dfz_dx) = self.xp1.map(|xp1| {
|
||||
// (xp1.y() - mid.y(), xp1.z() - mid.z())
|
||||
// }).unwrap_or_default();
|
||||
|
||||
// let (dfx_dy, dfz_dy) = self.yp1.map(|yp1| {
|
||||
// (yp1.x() - mid.x(), yp1.z() - mid.z())
|
||||
// }).unwrap_or_default();
|
||||
|
||||
// let (dfx_dz, dfy_dz) = self.zp1.map(|zp1| {
|
||||
// (zp1.x() - mid.x(), zp1.y() - mid.y())
|
||||
// }).unwrap_or_default();
|
||||
|
||||
let (dfy_dx, dfz_dx) = if self.xp1.is_some() {
|
||||
(self.xp1.unwrap().y() - mid.y(), self.xp1.unwrap().z() - mid.z())
|
||||
} else {
|
||||
(R::zero(), R::zero())
|
||||
};
|
||||
|
||||
let (dfx_dy, dfz_dy) = if self.yp1.is_some() {
|
||||
(self.yp1.unwrap().x() - mid.x(), self.yp1.unwrap().z() - mid.z())
|
||||
} else {
|
||||
(R::zero(), R::zero())
|
||||
};
|
||||
|
||||
let (dfx_dz, dfy_dz) = if self.zp1.is_some() {
|
||||
(self.zp1.unwrap().x() - mid.x(), self.zp1.unwrap().y() - mid.y())
|
||||
} else {
|
||||
(R::zero(), R::zero())
|
||||
};
|
||||
|
||||
FieldDeltas {
|
||||
dfy_dx,
|
||||
dfz_dx,
|
||||
dfx_dy,
|
||||
dfz_dy,
|
||||
dfx_dz,
|
||||
dfy_dz,
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
struct FieldDeltas<R> {
|
||||
dfy_dx: R,
|
||||
dfz_dx: R,
|
||||
dfx_dy: R,
|
||||
dfz_dy: R,
|
||||
dfx_dz: R,
|
||||
dfy_dz: R,
|
||||
}
|
||||
|
||||
impl<R: Real> FieldDeltas<R> {
|
||||
fn nabla(self) -> Vec3<R> {
|
||||
Vec3::new(
|
||||
self.dfz_dy - self.dfy_dz,
|
||||
self.dfx_dz - self.dfz_dx,
|
||||
self.dfy_dx - self.dfx_dy,
|
||||
)
|
||||
}
|
||||
}
|
||||
|
||||
pub struct StepEContext<'a, R, M> {
|
||||
pub inv_feature_size: R,
|
||||
pub time_step: R,
|
||||
pub stim_e: Vec3<R>,
|
||||
pub mat: &'a M,
|
||||
/// Input field sampled near this location
|
||||
pub in_h: VolumeSampleNeg<R>,
|
||||
pub in_e: Vec3<R>,
|
||||
}
|
||||
|
||||
impl<'a, R: Real, M: Material<R>> StepEContext<'a, R, M> {
|
||||
pub fn step_flat_view<RM, RF, WF>(
|
||||
meta: SimMeta<R>,
|
||||
mat: &RM,
|
||||
stim_e: &RF,
|
||||
e: &mut WF,
|
||||
h: &RF,
|
||||
idx: Vec3u,
|
||||
)
|
||||
where
|
||||
RM: Index<usize, Output=M> + ?Sized,
|
||||
RF: Index<usize, Output=Vec3<R>> + ?Sized,
|
||||
WF: Index<usize, Output=Vec3<R>> + IndexMut<usize> + ?Sized,
|
||||
{
|
||||
let dim = meta.dim;
|
||||
let stim_e_matrix = DimensionedSlice::new(dim, stim_e);
|
||||
let mat_matrix = DimensionedSlice::new(dim, mat);
|
||||
let mut e_matrix = DimensionedSlice::new(dim, e);
|
||||
let h_matrix = DimensionedSlice::new(dim, h);
|
||||
|
||||
let stim_e = stim_e_matrix[idx];
|
||||
let mat = &mat_matrix[idx];
|
||||
let in_e = e_matrix[idx];
|
||||
let in_h = VolumeSampleNeg::from_indexable(&h_matrix, idx);
|
||||
|
||||
let update_state = StepEContext {
|
||||
inv_feature_size: meta.inv_feature_size,
|
||||
time_step: meta.time_step,
|
||||
stim_e,
|
||||
mat,
|
||||
in_h,
|
||||
in_e,
|
||||
};
|
||||
let new_e = update_state.step_e();
|
||||
e_matrix[idx] = new_e;
|
||||
}
|
||||
|
||||
pub fn step_e(self) -> Vec3<R> {
|
||||
let twice_eps0 = R::twice_eps0();
|
||||
let deltas = self.in_h.delta_h();
|
||||
// \nabla x H
|
||||
let nabla_h = deltas.nabla() * self.inv_feature_size;
|
||||
// $\nabla x H = \epsilon_0 dE/dt + \sigma E$
|
||||
// no-conductivity version:
|
||||
// let delta_e = nabla_h * (self.time_step * EPS0_INV);
|
||||
let sigma = self.mat.conductivity();
|
||||
let e_prev = self.in_e;
|
||||
let delta_e = (nabla_h - e_prev.elem_mul(sigma)).elem_div(
|
||||
sigma*self.time_step + Vec3::uniform(twice_eps0)
|
||||
)*(R::two()*self.time_step);
|
||||
// println!("spirv-step_e delta_e: {:?}", delta_e);
|
||||
e_prev + delta_e + self.stim_e
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
pub struct StepHContext<'a, R, M> {
|
||||
pub inv_feature_size: R,
|
||||
pub time_step: R,
|
||||
pub stim_h: Vec3<R>,
|
||||
pub mat: &'a M,
|
||||
/// Input field sampled near this location
|
||||
pub in_e: VolumeSamplePos<R>,
|
||||
pub in_h: Vec3<R>,
|
||||
pub in_m: Vec3<R>,
|
||||
}
|
||||
|
||||
impl<'a, R: Real, M: Material<R>> StepHContext<'a, R, M> {
|
||||
pub fn step_flat_view<RM, RF, WF>(
|
||||
meta: SimMeta<R>,
|
||||
mat: &RM,
|
||||
stim_h: &RF,
|
||||
e: &RF,
|
||||
h: &mut WF,
|
||||
m: &mut WF,
|
||||
idx: Vec3u,
|
||||
)
|
||||
where
|
||||
RM: Index<usize, Output=M> + ?Sized,
|
||||
RF: Index<usize, Output=Vec3<R>> + ?Sized,
|
||||
WF: Index<usize, Output=Vec3<R>> + IndexMut<usize> + ?Sized,
|
||||
{
|
||||
let dim = meta.dim;
|
||||
let stim_h_matrix = DimensionedSlice::new(dim, stim_h);
|
||||
let mat_matrix = DimensionedSlice::new(dim, mat);
|
||||
let e_matrix = DimensionedSlice::new(dim, e);
|
||||
let mut h_matrix = DimensionedSlice::new(dim, h);
|
||||
let mut m_matrix = DimensionedSlice::new(dim, m);
|
||||
|
||||
let stim_h = stim_h_matrix[idx];
|
||||
let mat = &mat_matrix[idx];
|
||||
let in_e = VolumeSamplePos::from_indexable(&e_matrix, dim, idx);
|
||||
let in_h = h_matrix[idx];
|
||||
let in_m = m_matrix[idx];
|
||||
|
||||
let update_state = StepHContext {
|
||||
inv_feature_size: meta.inv_feature_size,
|
||||
time_step: meta.time_step,
|
||||
stim_h,
|
||||
mat,
|
||||
in_e,
|
||||
in_h,
|
||||
in_m,
|
||||
};
|
||||
let (new_h, new_m) = update_state.step_h();
|
||||
h_matrix[idx] = new_h;
|
||||
m_matrix[idx] = new_m;
|
||||
}
|
||||
|
||||
pub fn step_h(self) -> (Vec3<R>, Vec3<R>) {
|
||||
let mu0 = R::mu0();
|
||||
let mu0_inv = R::mu0_inv();
|
||||
let deltas = self.in_e.delta_e();
|
||||
// println!("spirv-step_h delta_e_struct: {:?}", deltas);
|
||||
// \nabla x E
|
||||
let nabla_e = deltas.nabla() * self.inv_feature_size;
|
||||
// println!("spirv-step_h nabla_e: {:?}", nabla_e);
|
||||
let delta_b = nabla_e * (-self.time_step);
|
||||
|
||||
// Relation between these is: B = mu0*(H + M)
|
||||
let old_h = self.in_h;
|
||||
let old_m = self.in_m;
|
||||
let old_b = (old_h + old_m) * mu0;
|
||||
|
||||
let new_b = old_b + delta_b + self.stim_h * mu0;
|
||||
let mat = self.mat;
|
||||
let new_m = mat.move_b_vec(old_m, new_b);
|
||||
let new_h = new_b * mu0_inv - new_m;
|
||||
// println!("spirv-step_h delta_h: {:?}", delta_h);
|
||||
(new_h, new_m)
|
||||
}
|
||||
}
|
||||
Reference in New Issue
Block a user