use crate::geom::{HasCrossSection, Meters, Region, Torus, WorldRegion}; use crate::real::{Real as _, ToFloat as _}; use crate::cross::vec::{Vec3, Vec3u}; use crate::sim::AbstractSim; use serde::{Serialize, Deserialize}; use std::ops::AddAssign; // TODO: do we really need both Send and Sync? pub trait AbstractMeasurement: Send + Sync { fn key_value(&self, state: &S) -> Vec; } pub fn as_dyn_measurements>(meas: &[M]) -> Vec<&dyn AbstractMeasurement> { meas.into_iter().map(|m| m as &dyn AbstractMeasurement).collect() } /// combine several measurements pub fn eval_multiple(state: &S, meas: &[&dyn AbstractMeasurement]) -> Vec { meas.into_iter().flat_map(|m| m.key_value(state).into_iter()).collect() } #[derive(Clone, Copy, Debug, PartialEq, Serialize, Deserialize)] pub enum MeasurementValue { Field(Vec3), Float(f32), Int(u64), Dim(Vec3u), } impl From> for MeasurementValue { fn from(v: Vec3) -> Self { Self::Field(v) } } impl From for MeasurementValue { fn from(v: f32) -> Self { Self::Float(v) } } impl From for MeasurementValue { fn from(v: u64) -> Self { Self::Int(v) } } impl From for MeasurementValue { fn from(v: Vec3u) -> Self { Self::Dim(v) } } #[derive(Clone, Debug, PartialEq, Serialize, Deserialize)] pub struct Measurement { name: String, value: MeasurementValue, /// e.g. "A" for Amps unit: String, } impl Measurement { fn new>(name: &str, value: T, unit: &str) -> Self { Self { name: name.to_owned(), value: value.into(), unit: unit.to_owned(), } } fn new_unitless>(name: &str, value: T) -> Self { Self::new(name, value, "") } pub fn name(&self) -> &str { &self.name } pub fn pretty_print(&self) -> String { use MeasurementValue::*; match self.value { Field(v) => format!("{}{}", v, self.unit), Float(f) => if self.unit != "" { SiScale::format_short(f, &self.unit) } else { f.to_string() }, Int(u) => format!("{}{}", u, self.unit), Dim(v) => format!("{}x{}x{}{}", v.x(), v.y(), v.z(), self.unit), } } /// format the Measurement in a way that could be parseable later. /// one major use case for this is in dumping the type to a CSV. pub fn machine_readable(&self) -> String { use MeasurementValue::*; match self.value { Field(v) => format!("{},{},{}", v.x(), v.y(), v.z()), Float(f) => f.to_string(), Int(u) => u.to_string(), Dim(v) => format!("{},{},{}", v.x(), v.y(), v.z()), } } /// retrieve the float value of this measurement -- if it's of float type. /// useful for tests pub fn get_float(&self) -> Option { match self.value { MeasurementValue::Float(f) => Some(f), _ => None, } } } impl AbstractMeasurement for Measurement { fn key_value(&self, _state: &S) -> Vec { vec![self.clone()] } } enum SiScale { Pico, Nano, Micro, Milli, Unit, Kilo, Mega, Giga, Terra, } impl SiScale { fn for_value(v: f32) -> Self { use SiScale::*; match v.abs() { v if v < 1e-12 => Unit, v if v < 1e-9 => Pico, v if v < 1e-6 => Nano, v if v < 1e-3 => Micro, v if v < 1e0 => Milli, v if v < 1e3 => Unit, v if v < 1e6 => Kilo, v if v < 1e9 => Mega, v if v < 1e12 => Giga, v if v < 1e15 => Terra, _ => Unit } } /// return the numerical scale of this prefix. /// e.g. `scale(&Pico) -> 1e-12 fn scale(&self) -> f32 { use SiScale::*; match *self { Pico => 1e-12, Nano => 1e-9, Micro => 1e-6, Milli => 1e-3, Unit => 1.0, Kilo => 1e3, Mega => 1e6, Giga => 1e9, Terra => 1e12, } } /// return the short string for this scale. /// e.g. `shortcode(Pico) -> "p"` fn shortcode(&self) -> &'static str { use SiScale::*; match *self { Pico => "p", Nano => "n", Micro => "u", Milli => "m", Unit => "", Kilo => "k", Mega => "M", Giga => "G", Terra => "T", } } /// format `v`, with the provided unit. /// e.g. `format_short(1234, "A") -> "1.23 kA" fn format_short(v: f32, unit: &str) -> String { let si = SiScale::for_value(v); let scaled = v / si.scale(); format!("{:.2} {}{}", scaled, si.shortcode(), unit) } } #[derive(Clone, Serialize, Deserialize)] pub struct Time; impl AbstractMeasurement for Time { fn key_value(&self, state: &S) -> Vec { vec![ Measurement::new_unitless("step", state.step_no()), Measurement::new("time", state.time(), "s"), ] } } #[derive(Clone, Serialize, Deserialize)] pub struct Meta; impl AbstractMeasurement for Meta { fn key_value(&self, state: &S) -> Vec { vec![ Measurement::new_unitless("dim", state.size().0), Measurement::new("feature_size", state.feature_size(), "m"), ] } } pub struct Volume { name: String, region: Box, } impl Volume { pub fn new(name: &str, r: R) -> Self { Self { name: name.into(), region: Box::new(r) } } /// Returns the volume of the region, in units of um^3 fn data(&self, state: &S) -> f32 { let feat_um = state.feature_size() as f64 * 1e6; (state.volume_of_region(&*self.region) as f64 * feat_um * feat_um * feat_um) as f32 } } impl AbstractMeasurement for Volume { fn key_value(&self, state: &S) -> Vec { vec![ Measurement::new(&format!("Vol({})", self.name), self.data(state), "um^3"), ] } } pub struct Current { name: String, region: Box, } impl Current { pub fn new(name: &str, r: R) -> Self { Self { name: name.into(), region: Box::new(r) } } fn data(&self, state: &S) -> (f32, Vec3) { let FieldSample(volume, current_mag, current_vec) = state.map_sum_over_enumerated(&*self.region, |coord: Meters, _cell| { let current = state.current(coord); FieldSample(1, current.mag().cast(), current.cast()) }); let mean_current_mag = current_mag.to_f32() / (f32::from_primitive(volume)); let mean_current_vec = current_vec.cast::() / (f32::from_primitive(volume)); (mean_current_mag, mean_current_vec.cast()) } } // TODO: clean up these FieldSample types #[derive(Default)] struct TupleSum(T); impl std::iter::Sum for TupleSum<(T0, T1)> { fn sum(iter: I) -> Self where I: Iterator { let mut s = Self::default(); for TupleSum((a0, a1)) in iter { s.0.0 += a0; s.0.1 += a1; } s } } #[derive(Default)] struct FieldSample(u32, f64, Vec3); impl std::iter::Sum for FieldSample { fn sum(iter: I) -> Self where I: Iterator { let mut s = FieldSample::default(); for FieldSample(a, b, c) in iter { s = FieldSample(s.0 + a, s.1 + b, s.2 + c); } s } } #[derive(Default)] struct FieldSamples(T); impl std::iter::Sum for FieldSamples<[FieldSample; 2]> { fn sum(iter: I) -> Self where I: Iterator { let mut s = Self::default(); for p in iter { s.0[0] = FieldSample(s.0[0].0 + p.0[0].0, s.0[0].1 + p.0[0].1, s.0[0].2 + p.0[0].2); s.0[1] = FieldSample(s.0[1].0 + p.0[1].0, s.0[1].1 + p.0[1].1, s.0[1].2 + p.0[1].2); } s } } impl std::iter::Sum for FieldSamples<[FieldSample; 3]> { fn sum(iter: I) -> Self where I: Iterator { let mut s = Self::default(); for p in iter { s.0[0] = FieldSample(s.0[0].0 + p.0[0].0, s.0[0].1 + p.0[0].1, s.0[0].2 + p.0[0].2); s.0[1] = FieldSample(s.0[1].0 + p.0[1].0, s.0[1].1 + p.0[1].1, s.0[1].2 + p.0[1].2); s.0[2] = FieldSample(s.0[2].0 + p.0[2].0, s.0[2].1 + p.0[2].1, s.0[2].2 + p.0[2].2); } s } } impl AbstractMeasurement for Current { fn key_value(&self, state: &S) -> Vec { let (mean_current_mag, mean_current_vec) = self.data(state); vec![ Measurement::new( &format!("Imag/cell({})", self.name), mean_current_mag, "A", ), Measurement::new( &format!("I/cell({})", self.name), mean_current_vec, "A", ), ] } } /// Measures the current directed around a closed loop #[derive(Clone, Serialize, Deserialize)] pub struct CurrentLoop { name: String, region: R, } impl CurrentLoop { pub fn new(name: &str, r: R) -> Self { Self { name: name.into(), region: r, } } } impl CurrentLoop { fn data(&self, state: &S) -> f32 { // - current exists as a property of a 2d surface. // - the user has provided us a 3d volume which behaves as though it's an extruded surface: // for any point in the volume we can query the normal vector of the cross section // containing that point. // - we choose that measuring the "current" on such a volume means to measure the average // current through all its cross sections (for most boring materials, each // cross section has nearly identical current). // - therefore, enumerate the entire volume and compute the "net" current (the sum over // each cell of whatever current in that cell is along the cross-section normal). // then divide by the number of complete cross sections we measured, to average. let feature_area = state.feature_size() * state.feature_size(); let TupleSum((net_current, cross_sections)) = state.map_sum_over_enumerated(&self.region, move |coord: Meters, _cell| { // `normal` represents both the size of the cross section (m^2) this cell belongs to, // and the normal direction of the cross section. let normal = self.region.cross_section_normal(coord); // [m^2] // calculate the amount of normal current through this specific cell let current_density = state.current_density(coord); // [A/m^2] let cross_sectional_current = feature_area * current_density.dot(normal.norm()); // [A] // keep track of how many cross sections we enumerate, since each additional cross // sections represents a double-count of the current. let num_cross_sections_filled = feature_area / normal.mag(); TupleSum((cross_sectional_current, num_cross_sections_filled)) }); let mean_cross_sectional_current = net_current.cast::() / cross_sections; mean_cross_sectional_current } } impl AbstractMeasurement for CurrentLoop { fn key_value(&self, state: &S) -> Vec { let cross_sectional_current = self.data(state); vec![ Measurement::new( &format!("I({})", self.name), cross_sectional_current, "A" ), ] } } /// Measures the M, B field directed around a closed loop #[derive(Clone, Serialize, Deserialize)] pub struct MagneticLoop { name: String, region: Torus } impl MagneticLoop { pub fn new(name: &str, r: Torus) -> Self { Self { name: name.into(), region: r, } } fn data(&self, state: &S) -> (f32, f32, f32) { let FieldSamples([ FieldSample(volume, directed_m, _m_vec), FieldSample(_, directed_b, _b_vec), FieldSample(_, directed_h, _h_vec), ]) = state.map_sum_over_enumerated(&self.region, |coord: Meters, cell| { let normal = self.region.axis(); let to_coord = *coord - *self.region.center(); let tangent = normal.cross(to_coord).norm(); let m = cell.m(); let directed_m = m.dot(tangent.cast()); let b = cell.b(); let directed_b = b.dot(tangent.cast()); let h = cell.h(); let directed_h = h.dot(tangent.cast()); FieldSamples([ FieldSample(1, directed_m.cast(), m.cast()), FieldSample(1, directed_b.cast(), b.cast()), FieldSample(1, directed_h.cast(), h.cast()), ]) }); // let cross_section = self.region.cross_section() / (state.feature_size() * state.feature_size()); let mean_directed_m = directed_m.cast::() / f32::from_primitive(volume); // let cross_sectional_m = mean_directed_m * cross_section; let mean_directed_b = directed_b.cast::() / f32::from_primitive(volume); // let cross_sectional_b = mean_directed_b * cross_section; let mean_directed_h = directed_h.cast::() / f32::from_primitive(volume); // let cross_sectional_h = mean_directed_h * cross_section; // format!( // "M({}): {:.2e}; B({}): {:.2e}; H({}): {:.2e}", // self.name, cross_sectional_m, // self.name, cross_sectional_b, // self.name, cross_sectional_h, // ) (mean_directed_m, mean_directed_b, mean_directed_h) } } impl AbstractMeasurement for MagneticLoop { fn key_value(&self, state: &S) -> Vec { let (mean_directed_m, mean_directed_b, mean_directed_h) = self.data(state); vec![ Measurement::new_unitless(&format!("M({})", self.name), mean_directed_m), Measurement::new_unitless(&format!("B({})", self.name), mean_directed_b), Measurement::new_unitless(&format!("H({})", self.name), mean_directed_h), ] } } /// mean M over a region pub struct MagneticFlux { name: String, region: Box, } impl MagneticFlux { pub fn new(name: &str, r: R) -> Self { Self { name: name.into(), region: Box::new(r) } } fn data(&self, state: &S) -> Vec3 { let FieldSample(volume, _directed_mag, mag_vec) = state.map_sum_over(&*self.region, |cell| { let b = cell.b(); let mag = b.mag(); FieldSample(1, mag.cast(), b.cast()) }); let mean_mag = mag_vec.cast() / f32::from_primitive(volume); mean_mag } } impl AbstractMeasurement for MagneticFlux { fn key_value(&self, state: &S) -> Vec { let mean_mag = self.data(state); vec![ Measurement::new_unitless( &format!("Bavg({})", self.name), mean_mag, ) ] } } /// mean B over a region pub struct Magnetization { name: String, region: Box, } impl Magnetization { pub fn new(name: &str, r: R) -> Self { Self { name: name.into(), region: Box::new(r) } } fn data(&self, state: &S) -> Vec3 { let FieldSample(volume, _directed_mag, mag_vec) = state.map_sum_over(&*self.region, |cell| { let m = cell.m(); let mag = m.mag(); FieldSample(1, mag.cast(), m.cast()) }); let mean_mag = mag_vec.cast() / f32::from_primitive(volume); mean_mag } } impl AbstractMeasurement for Magnetization { fn key_value(&self, state: &S) -> Vec { let mean_mag = self.data(state); vec![ Measurement::new_unitless( &format!("Mavg({})", self.name), mean_mag ), ] } } fn loc(v: Meters) -> String { format!("{:.0} um", *v * f32::from_primitive(1_000_000)) } /// M #[derive(Clone, Serialize, Deserialize)] pub struct MagnetizationAt(pub Meters); impl AbstractMeasurement for MagnetizationAt { fn key_value(&self, state: &S) -> Vec { let m = state.sample(self.0).m(); vec![ Measurement::new_unitless(&format!("M{}", loc(self.0)), m.cast()) ] } } /// B #[derive(Clone, Serialize, Deserialize)] pub struct MagneticFluxAt(pub Meters); impl AbstractMeasurement for MagneticFluxAt { fn key_value(&self, state: &S) -> Vec { let b = state.sample(self.0).b(); vec![ Measurement::new_unitless( &format!("B{}", loc(self.0)), b.cast() ) ] } } /// H #[derive(Clone, Serialize, Deserialize)] pub struct MagneticStrengthAt(pub Meters); impl AbstractMeasurement for MagneticStrengthAt { fn key_value(&self, state: &S) -> Vec { let h = state.sample(self.0).h(); vec![ Measurement::new_unitless( &format!("H{}", loc(self.0)), h.cast() ) ] } } #[derive(Clone, Serialize, Deserialize)] pub struct ElectricField(pub Meters); impl AbstractMeasurement for ElectricField { fn key_value(&self, state: &S) -> Vec { let e = state.sample(self.0).e(); vec![ Measurement::new_unitless( &format!("E{}", loc(self.0)), e.cast() ) ] } } pub struct Energy { name: String, region: Box, } impl Energy { pub fn world() -> Self { Self::new("World", WorldRegion) } pub fn new(name: &str, region: R) -> Self { Self { name: name.into(), region: Box::new(region), } } pub(crate) fn data(&self, state: &S) -> f32 { // Potential energy stored in a E/M field: // https://en.wikipedia.org/wiki/Magnetic_energy // https://en.wikipedia.org/wiki/Electric_potential_energy#Energy_stored_in_an_electrostatic_field_distribution // TODO: consider the M field? https://en.wikipedia.org/wiki/Potential_energy#Magnetic_potential_energy // U(B) = 1/2 \int H . B dV // U(E) = 1/2 \int E . D dV #[allow(non_snake_case)] let dV = state.feature_volume(); let e = f64::from_primitive(0.5 * dV) * state.map_sum_over(&*self.region, |cell| { // E . D = E . (E + P) = E.E since we don't model polarization fields cell.h().dot(cell.b()).to_f64() + cell.e().mag_sq().to_f64() }); e.cast() } } impl AbstractMeasurement for Energy { fn key_value(&self, state: &S) -> Vec { let e = self.data(state); vec![ Measurement::new( &format!("U({})", self.name), e, "J" ) ] } } pub struct Power { name: String, region: Box } impl Power { pub fn world() -> Self { Self::new("World", WorldRegion) } pub fn new(name: &str, region: R) -> Self { Self { name: name.into(), region: Box::new(region), } } fn data(&self, state: &S) -> f32 { // Power is P = IV = A*J*V = L^2*J.(LE) = L^3 J.E // where L is feature size. #[allow(non_snake_case)] let dV = state.feature_volume(); let power = f64::from_primitive(dV) * state.map_sum_over(&*self.region, |cell| { cell.current_density().dot(cell.e()).to_f64() }); power.cast() } } impl AbstractMeasurement for Power { fn key_value(&self, state: &S) -> Vec { let power = self.data(state); vec![ Measurement::new( &format!("P({})", self.name), power, "W" ) ] } } #[cfg(test)] pub mod test { use super::*; use crate::cross::mat::AnisomorphicConductor; use crate::cross::step::SimMeta; use crate::geom::Index; use crate::sim::{Fields, GenericSim}; use crate::stim::Stimulus; struct MockSim { e_field: Vec3, dim: Vec3u, feature_size: f32, mat: AnisomorphicConductor, } impl AbstractSim for MockSim { type Real = f32; type Material = AnisomorphicConductor; fn meta(&self) -> SimMeta { SimMeta::new(self.dim, self.feature_size, 1e-9) } fn step_no(&self) -> u64 { unimplemented!() } fn fields_at_index(&self, _pos: Index) -> Fields { Fields::new(self.e_field, Vec3::zero(), Vec3::zero()) } fn get_material_index(&self, _at: Index) -> &Self::Material { &self.mat } fn put_material_index(&mut self, _at: Index, _m: Self::Material) { unimplemented!() } fn step_multiple>(&mut self, _num_steps: u32, _s: &S) { unimplemented!() } fn to_generic(&self) -> GenericSim { unimplemented!() } } struct MockRegion { normal: Vec3, } impl HasCrossSection for MockRegion { fn cross_section_normal(&self, _p: Meters) -> Vec3 { self.normal } } impl Region for MockRegion { fn contains(&self, _p: Meters) -> bool { true } } #[test] fn current_loop_trivial() { let sim = MockSim { e_field: Vec3::new(1.0, 0.0, 0.0), dim: Vec3u::new(1, 1, 1), feature_size: 1.0, mat: AnisomorphicConductor::new(Vec3::new(1.0, 1.0, 1.0)), }; let region = MockRegion { normal: Vec3::new(1.0, 0.0, 0.0), }; let kv = CurrentLoop::new("test", region).key_value(&sim); assert_eq!(kv.len(), 1); // measured area is 1 m^2 // region cross-section is 1 m^2 // conductivity is 1 S/m assert_eq!(kv[0].get_float().unwrap(), 1.0); } #[test] fn current_loop_multi_cell() { let sim = MockSim { e_field: Vec3::new(1.0, 0.0, 0.0), dim: Vec3u::new(4, 4, 4), feature_size: 0.25, mat: AnisomorphicConductor::new(Vec3::new(1.0, 1.0, 1.0)), }; let region = MockRegion { normal: Vec3::new(1.0, 0.0, 0.0), }; let kv = CurrentLoop::new("test", region).key_value(&sim); assert_eq!(kv.len(), 1); // measured area is 1 m^2 // region cross-section is 1 m^2 // conductivity is 1 S/m assert_eq!(kv[0].get_float().unwrap(), 1.0); } #[test] fn current_loop_off_conductor() { let sim = MockSim { e_field: Vec3::new(1.0, 1.0, 1.0), dim: Vec3u::new(4, 4, 4), feature_size: 0.25, mat: AnisomorphicConductor::new(Vec3::new(0.0, 1.0, 1.0)), }; let region = MockRegion { normal: Vec3::new(1.0, 0.0, 0.0), }; let kv = CurrentLoop::new("test", region).key_value(&sim); assert_eq!(kv.len(), 1); // material is not conductive in the direction being queried assert_eq!(kv[0].get_float().unwrap(), 0.0); } #[test] fn current_loop_e_field() { let sim = MockSim { e_field: Vec3::new(4.0, 2.0, 1.0), dim: Vec3u::new(4, 4, 4), feature_size: 0.25, mat: AnisomorphicConductor::new(Vec3::new(1.0, 1.0, 1.0)), }; let region = MockRegion { normal: Vec3::new(1.0, 0.0, 0.0), }; let kv = CurrentLoop::new("test", region).key_value(&sim); assert_eq!(kv.len(), 1); // measured area is 1 m^2 // region cross-section is 1 m^2 // conductivity is 1 S/m // e field is 4 V/m assert_eq!(kv[0].get_float().unwrap(), 4.0); } #[test] fn current_loop_conductivity() { let sim = MockSim { e_field: Vec3::new(4.0, 2.0, 1.0), dim: Vec3u::new(4, 4, 4), feature_size: 0.25, mat: AnisomorphicConductor::new(Vec3::new(3.0, 1.0, 1.0)), }; let region = MockRegion { normal: Vec3::new(1.0, 0.0, 0.0), }; let kv = CurrentLoop::new("test", region).key_value(&sim); assert_eq!(kv.len(), 1); // measured area is 1 m^2 // region cross-section is 1 m^2 // conductivity is 3 S/m // e field is 4 V/m assert_eq!(kv[0].get_float().unwrap(), 3.0*4.0); } #[test] fn current_loop_cross_section() { let sim = MockSim { e_field: Vec3::new(4.0, 2.0, 1.0), dim: Vec3u::new(4, 4, 4), feature_size: 0.5, mat: AnisomorphicConductor::new(Vec3::new(3.0, 1.0, 1.0)), }; let region = MockRegion { normal: Vec3::new(16.0, 0.0, 0.0), }; let kv = CurrentLoop::new("test", region).key_value(&sim); assert_eq!(kv.len(), 1); // measured area is 2 m^2 // region cross-section is 16 m^2 // conductivity is 3 S/m // e field is 4 V/m assert_eq!(kv[0].get_float().unwrap(), 3.0*4.0*16.0); } }