diff --git a/crates/coremem/src/meas.rs b/crates/coremem/src/meas.rs index 83fccfb..01be8d4 100644 --- a/crates/coremem/src/meas.rs +++ b/crates/coremem/src/meas.rs @@ -3,6 +3,7 @@ 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 { @@ -243,6 +244,24 @@ impl Current { } } +// 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); @@ -323,33 +342,30 @@ impl CurrentLoop { } impl CurrentLoop { fn data(&self, state: &S) -> f32 { - // i use a statistical lens for this: - // 1. current is the rate of flow of charge into a surface. - // 2. in any context where it makes sense to think of current, the current through each - // cross-sectional **is the same**. - // 3. each point in our 3d region belongs to exactly one cross-sectional surface. - // 4. so, given a point: what's the expected current through the cross section it belongs to? - // - answer: that point's current density times the cross section's area. - // 5. average the above over the whole volume, and you get an "average current". - // - // we're sampling uniformly over the cell space -- not the set of cross sections. - // - however, if all cross sections have equal area, this is equivalent. - // sampling all points (instead of just a single point): - // 1) removes bias from step #4: current *within* a cross section is not uniform, but if - // we sample every point within the cross section and weight them equally, then the - // average is the truth. - // 2) probably combats grid quantization / artifacting. - let FieldSample(num_samples, sum_cross_sectional_current, _current_vec) = state.map_sum_over_enumerated(&self.region, |coord: Meters, _cell| { + // - 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] - // now we have an estimation of the entire current flowing through the cross section - // this cell belongs to. - let cross_sectional_current = current_density.dot(normal.cast()); // [A] - FieldSample(1, cross_sectional_current.cast(), current_density.cast()) + 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 = sum_cross_sectional_current.cast::() / f32::from_primitive(num_samples); + let mean_cross_sectional_current = net_current.cast::() / cross_sections; mean_cross_sectional_current } } diff --git a/crates/cross/src/real.rs b/crates/cross/src/real.rs index 2e2ab0c..62ad7e0 100644 --- a/crates/cross/src/real.rs +++ b/crates/cross/src/real.rs @@ -107,6 +107,10 @@ pub trait Real: self == Self::zero() } + fn inv(self) -> Self { + Self::one() / self + } + fn zero() -> Self; fn one() -> Self; fn two() -> Self;