Adom Fields

@@ -0,0 +1,195 @@+//! Partial-element equivalent circuit (PEEC) loop inductance from the board's copper (0.10.2, used by the pdn module).+//! The copper of two nets in a window is cut into cells on every copper layer; neighbouring cells of one net are+//! bars (x- or y-directed branches) and vias are z-directed branches between layers. Each branch has a partial self+//! inductance and a partial mutual inductance with every parallel branch (Ruehli: parallel filaments, the self at the+//! bar's geometric mean distance). At the frequencies where mounting inductance matters, jωL dwarfs the resistance,+//! so the inductive network is solved alone: [L Aᵀ; A 0][I; −φ] = [0; J], 1 A in at the port; the port's inductance+//! is its potential difference. A part's two pins can be shorted (the chip's rail pin to its ground pin, where the+//! decoupling loop closes inside the package). Validated against the 2-D field solve of a strip over a plane (tests).+use crate::board::Board;++const MU0_4PI: f64 = 1e-7; // H/m++/// mutual inductance of two parallel filaments along one axis: [a1, a2] and [b1, b2] (m), d apart (m)+pub fn m_par(a1: f64, a2: f64, b1: f64, b2: f64, d: f64) -> f64 {+    let f = |z: f64| z * (z / d).asinh() - (z * z + d * d).sqrt();+    MU0_4PI * (f(a2 - b1) - f(a2 - b2) - f(a1 - b1) + f(a1 - b2))+}+/// partial self inductance of a bar of length l, width w, thickness t (m): a filament at the cross-section's GMD+pub fn l_self(l: f64, w: f64, t: f64) -> f64 { m_par(0.0, l, 0.0, l, 0.2235 * (w + t)) }++#[derive(Clone, Copy, PartialEq)] pub enum Dir { X, Y, Z }+#[derive(Clone)] pub struct Branch { pub n1: usize, pub n2: usize, pub dir: Dir, pub a: (f64, f64, f64), pub len: f64, pub w: f64, pub t: f64 }++/// the loop inductance seen at a port (node p, node q) of a branch network: dense LU on the MNA system+pub fn solve(_nn: usize, br: &[Branch], shorts: &[(usize, usize)], p: usize, q: usize) -> Option<f64> {+    // only the nodes a branch or a short touches are unknowns (empty cells would leave the system singular)+    // keep only the piece of copper the port is on: floating islands leave their potentials undefined+    let mut adj: std::collections::BTreeMap<usize, Vec<usize>> = Default::default();+    for (a, b) in br.iter().map(|b| (b.n1, b.n2)).chain(shorts.iter().copied()) { adj.entry(a).or_default().push(b); adj.entry(b).or_default().push(a); }+    let mut keep: std::collections::BTreeSet<usize> = Default::default(); let mut stack = vec![p];+    while let Some(u) = stack.pop() { if keep.insert(u) { for v in adj.get(&u).into_iter().flatten() { if !keep.contains(v) { stack.push(*v); } } } }+    if !keep.contains(&q) { return None; }+    let br: Vec<Branch> = br.iter().filter(|b| keep.contains(&b.n1)).cloned().collect();+    let shorts: Vec<(usize, usize)> = shorts.iter().filter(|s| keep.contains(&s.0)).copied().collect();+    let mut used: std::collections::BTreeMap<usize, usize> = Default::default();+    for (a, b) in br.iter().map(|b| (b.n1, b.n2)).chain(shorts.iter().copied()) { let k = used.len(); used.entry(a).or_insert(k); let k = used.len(); used.entry(b).or_insert(k); }+    let (Some(&p), Some(&q)) = (used.get(&p), used.get(&q)) else { return None };+    let br: Vec<Branch> = br.iter().map(|b| Branch { n1: used[&b.n1], n2: used[&b.n2], ..b.clone() }).collect();+    let shorts: Vec<(usize, usize)> = shorts.iter().map(|(a, b)| (used[a], used[b])).collect();+    let (br, shorts, nn) = (&br[..], &shorts[..], used.len());+    let nb = br.len(); let ns = shorts.len();+    // unknowns: branch currents (nb), short currents (ns), node potentials except node q (ground) (nn - 1)+    let idx = |n: usize| -> Option<usize> { if n == q { None } else if n < q { Some(nb + ns + n) } else { Some(nb + ns + n - 1) } };+    let m = nb + ns + nn - 1;+    let mut a = vec![0.0f64; m * m]; let mut rhs = vec![0.0f64; m];+    for i in 0..nb { for j in i..nb {+        let (bi, bj) = (&br[i], &br[j]); if bi.dir != bj.dir { continue; }+        let v = if i == j { l_self(bi.len, bi.w, bi.t) } else {+            // the axis coordinate and the perpendicular distance of two parallel bars+            let (s1, s2, perp) = match bi.dir { Dir::X => (bi.a.0, bj.a.0, ((bi.a.1 - bj.a.1).powi(2) + (bi.a.2 - bj.a.2).powi(2)).sqrt()), Dir::Y => (bi.a.1, bj.a.1, ((bi.a.0 - bj.a.0).powi(2) + (bi.a.2 - bj.a.2).powi(2)).sqrt()), Dir::Z => (bi.a.2, bj.a.2, ((bi.a.0 - bj.a.0).powi(2) + (bi.a.1 - bj.a.1).powi(2)).sqrt()) };+            let d = perp.max(0.2235 * (bi.w + bi.t));+            m_par(s1, s1 + bi.len, s2, s2 + bj.len, d)+        };+        a[i * m + j] = v; a[j * m + i] = v;+    } }+    // KVL rows: L·I − (φ(n1) − φ(n2)) = 0 → L·I + φ(n2) − φ(n1); signs folded into the incidence+    let inc = |a: &mut Vec<f64>, row: usize, n1: usize, n2: usize| { if let Some(k) = idx(n1) { a[row * m + k] -= 1.0; } if let Some(k) = idx(n2) { a[row * m + k] += 1.0; } };+    for (i, b) in br.iter().enumerate() { inc(&mut a, i, b.n1, b.n2); }+    for (s, (n1, n2)) in shorts.iter().enumerate() { inc(&mut a, nb + s, *n1, *n2); }+    // KCL rows (one per node except q): currents leaving through branches and shorts = injected+    for (i, b) in br.iter().enumerate() { if let Some(r) = idx(b.n1) { a[r * m + i] += 1.0; } if let Some(r) = idx(b.n2) { a[r * m + i] -= 1.0; } }+    for (s, (n1, n2)) in shorts.iter().enumerate() { let c = nb + s; if let Some(r) = idx(*n1) { a[r * m + c] += 1.0; } if let Some(r) = idx(*n2) { a[r * m + c] -= 1.0; } }+    let rp = idx(p)?; rhs[rp] = 1.0; // 1 A into p, out at q+    let x = lu_solve(&mut a, &mut rhs, m)?;+    Some(x[rp])+}++fn lu_solve(a: &mut [f64], b: &mut [f64], m: usize) -> Option<Vec<f64>> {+    for k in 0..m {+        let piv = (k..m).max_by(|i, j| a[i * m + k].abs().partial_cmp(&a[j * m + k].abs()).unwrap())?;+        if a[piv * m + k].abs() < 1e-30 { return None; }+        if piv != k { for c in 0..m { a.swap(k * m + c, piv * m + c); } b.swap(k, piv); }+        let d = a[k * m + k];+        let (top, rest) = a.split_at_mut((k + 1) * m); let rowk = &top[k * m..(k + 1) * m];+        for i in 0..(m - k - 1) {+            let row = &mut rest[i * m..(i + 1) * m]; let f = row[k] / d; if f == 0.0 { continue; }+            row[k] = 0.0; for c in (k + 1)..m { row[c] -= f * rowk[c]; } b[k + 1 + i] -= f * b[k];+        }+    }+    let mut x = vec![0.0; m];+    for k in (0..m).rev() { let mut s = b[k]; for c in (k + 1)..m { s -= a[k * m + c] * x[c]; } x[k] = s / a[k * m + k]; }+    Some(x)+}++/// A grid of copper cells: for each layer and cell, the net there (or none). Builds the branch network.+pub struct Grid { pub x0: f64, pub y0: f64, pub d: f64, pub w: usize, pub h: usize, pub z: Vec<f64>, pub t: f64, pub cell: Vec<Vec<Option<u16>>> }+impl Grid {+    pub fn node(&self, k: usize, r: usize, c: usize) -> usize { (k * self.h + r) * self.w + c }+    /// branches between neighbouring cells of one net on one layer, and z branches at the given via points (layer pairs)+    pub fn network(&self, vias: &[(f64, f64, u16)]) -> Vec<Branch> {+        let mut br = Vec::new(); let (d, t) = (self.d * 1e-3, self.t);+        for k in 0..self.cell.len() { for r in 0..self.h { for c in 0..self.w {+            let Some(n) = self.cell[k][r * self.w + c] else { continue };+            let (x, y, z) = ((self.x0 + (c as f64 + 0.5) * self.d) * 1e-3, (self.y0 + (r as f64 + 0.5) * self.d) * 1e-3, self.z[k] * 1e-3);+            if c + 1 < self.w && self.cell[k][r * self.w + c + 1] == Some(n) { br.push(Branch { n1: self.node(k, r, c), n2: self.node(k, r, c + 1), dir: Dir::X, a: (x, y, z), len: d, w: d, t }); }+            if r + 1 < self.h && self.cell[k][(r + 1) * self.w + c] == Some(n) { br.push(Branch { n1: self.node(k, r, c), n2: self.node(k, r + 1, c), dir: Dir::Y, a: (x, y, z), len: d, w: d, t }); }+        } } }+        for (vx, vy, n) in vias {+            let (c, r) = (((vx - self.x0) / self.d) as isize, ((vy - self.y0) / self.d) as isize); if c < 0 || r < 0 || c as usize >= self.w || r as usize >= self.h { continue; }+            let (r, c) = (r as usize, c as usize);+            let ks: Vec<usize> = (0..self.cell.len()).filter(|k| self.cell[*k][r * self.w + c] == Some(*n)).collect();+            for p in ks.windows(2) { let (za, zb) = (self.z[p[0]], self.z[p[1]]); br.push(Branch { n1: self.node(p[0], r, c), n2: self.node(p[1], r, c), dir: Dir::Z, a: ((self.x0 + (c as f64 + 0.5) * self.d) * 1e-3, (self.y0 + (r as f64 + 0.5) * self.d) * 1e-3, za.min(zb) * 1e-3), len: (za - zb).abs() * 1e-3, w: 0.3e-3, t: 0.3e-3 }); }+        }+        br+    }+}++/// Rasterise the board's copper for the named nets inside a window (mm) at pitch d (mm).+pub fn grid_from_board(b: &Board, nets: &[String], win: [f64; 4], d: f64, z_layers: &[(String, f64)], t_mm: f64) -> Grid {+    let w = ((win[2] - win[0]) / d).ceil() as usize; let h = ((win[3] - win[1]) / d).ceil() as usize;+    let inside = |p: (f64, f64), poly: &[(f64, f64)]| { let mut c = false; let n = poly.len(); for i in 0..n { let (a, bb) = (poly[i], poly[(i + n - 1) % n]); if (a.1 > p.1) != (bb.1 > p.1) && p.0 < (bb.0 - a.0) * (p.1 - a.1) / (bb.1 - a.1) + a.0 { c = !c; } } c };+    let seg = |p: (f64, f64), a: (f64, f64), bb: (f64, f64)| { let (dx, dy) = (bb.0 - a.0, bb.1 - a.1); let l2 = dx * dx + dy * dy; let t = if l2 == 0.0 { 0.0 } else { (((p.0 - a.0) * dx + (p.1 - a.1) * dy) / l2).clamp(0.0, 1.0) }; (p.0 - a.0 - t * dx).hypot(p.1 - a.1 - t * dy) };+    let mut cell: Vec<Vec<Option<u16>>> = vec![vec![None; w * h]; z_layers.len()];+    for (k, (ln, _)) in z_layers.iter().enumerate() { for r in 0..h { for c in 0..w {+        let p = (win[0] + (c as f64 + 0.5) * d, win[1] + (r as f64 + 0.5) * d);+        for (ni, n) in nets.iter().enumerate() {+            let hit = b.polys.iter().chain(b.pad_polys.iter()).any(|q| q.net == *n && (q.layer == *ln || q.layer == "*.Cu") && q.pts.len() > 2 && inside(p, &q.pts))+                || b.tracks.iter().any(|t| t.net == *n && t.layer == *ln && seg(p, t.a, t.b) <= t.w / 2.0);+            if hit { cell[k][r * w + c] = Some(ni as u16); break; }+        }+    } } }+    // copper narrower than a cell must not vanish (a 0.25 mm trace at 0.6 mm cells would break the loop): every cell a+    // track passes through, and the cell under every pad and via, belongs to its net unless another net holds it+    let mut mark = |k: usize, x: f64, y: f64, ni: u16, cell: &mut Vec<Vec<Option<u16>>>| { let (c, r) = (((x - win[0]) / d).floor(), ((y - win[1]) / d).floor()); if c >= 0.0 && r >= 0.0 && (c as usize) < w && (r as usize) < h { let i = r as usize * w + c as usize; if cell[k][i].is_none() { cell[k][i] = Some(ni); } } };+    for (k, (ln, _)) in z_layers.iter().enumerate() { for (ni, n) in nets.iter().enumerate() {+        for t in b.tracks.iter().filter(|t| t.net == *n && t.layer == *ln) { let len = (t.b.0 - t.a.0).hypot(t.b.1 - t.a.1); let steps = (len / (d / 3.0)).ceil().max(1.0) as usize; for i in 0..=steps { let f = i as f64 / steps as f64; mark(k, t.a.0 + f * (t.b.0 - t.a.0), t.a.1 + f * (t.b.1 - t.a.1), ni as u16, &mut cell); } }+        for p in b.pads.iter().filter(|p| p.net == *n && (p.layers.iter().any(|l| l == ln || l == "*.Cu"))) { mark(k, p.x, p.y, ni as u16, &mut cell); }+        for v in b.vias.iter().chain(b.plated_holes.iter()).filter(|v| v.net == *n) { mark(k, v.x, v.y, ni as u16, &mut cell); }+    } }+    Grid { x0: win[0], y0: win[1], d, w, h, z: z_layers.iter().map(|x| x.1).collect(), t: t_mm * 1e-3, cell }+}++/// The loop inductance (nH) between two nets: 1 A in at port_a (net a), out at port_b (net b), with a short from+/// short_a (net a) to short_b (net b) closing the loop, over the copper in the window. None when the loop is not+/// closed by the copper in the window, or the window would take more than max_branches.+pub fn loop_nh(b: &Board, layers: &[(String, f64)], t_mm: f64, nets: [&str; 2], port: [(f64, f64); 2], short: [(f64, f64); 2], win: [f64; 4], pitch: f64, max_branches: usize, resolve_mm: f64) -> Result<(f64, usize, f64), String> {+    let nets_v = vec![nets[0].to_string(), nets[1].to_string()];+    let vias: Vec<(f64, f64, u16)> = b.vias.iter().chain(b.plated_holes.iter()).filter_map(|v| Some((v.x, v.y, nets_v.iter().position(|q| *q == v.net)? as u16))).collect();+    // the dense solve grows as the cube of the branches: coarsen the cells until the window fits the budget (the+    // rasterising is cheap; at 0.5 mm the buck's loops are within about 5 % of their converged values)+    let (mut pitch, mut tries) = (pitch, 0);+    let (g, br) = loop {+        let g = grid_from_board(b, &nets_v, win, pitch, layers, t_mm); let br = g.network(&vias);+        if br.len() <= max_branches || tries >= 4 { break (g, br); }+        pitch *= ((br.len() as f64 / max_branches as f64).sqrt() * 1.05).max(1.1); tries += 1;+    };+    if br.is_empty() { return Err("no copper of the two nets in the window".into()); }+    // a plane closer than the cells can resolve gives a wrong loop (the validation: −37 % at cells 2.5 times the dielectric, −17 % at+    // 1.25, a few % at 1 or under): a thin four-layer stackup needs a faster solver than this dense one+    if pitch > resolve_mm { return Err(format!("the nearest plane is {resolve_mm:.2} mm away, closer than the {pitch:.2} mm cells this solve can afford in the window can resolve (thin stackups need a faster solver, not yet built)")); }+    if br.len() > max_branches || pitch > 1.0 { return Err(format!("the window needs {} branches even at {pitch:.2} mm cells", br.len())); }+    // a point's node: the cell under it on the first layer (top down) where that net's copper is+    let at = |p: (f64, f64), ni: u16| -> Option<usize> { let (c, r) = (((p.0 - g.x0) / g.d) as isize, ((p.1 - g.y0) / g.d) as isize); if c < 0 || r < 0 || c as usize >= g.w || r as usize >= g.h { return None; } let (r, c) = (r as usize, c as usize); (0..g.cell.len()).find(|k| g.cell[*k][r * g.w + c] == Some(ni)).map(|k| g.node(k, r, c)) };+    let (Some(p), Some(q)) = (at(port[0], 0), at(port[1], 1)) else { return Err("a capacitor pad is not on its net's copper in the grid".into()) };+    let (Some(sa), Some(sz)) = (at(short[0], 0), at(short[1], 1)) else { return Err("a load or ground pin is not on its net's copper in the grid".into()) };+    let l = solve(2 * g.w * g.h * layers.len(), &br, &[(sa, sz)], p, q).ok_or("the loop does not close through the copper in the window")?;+    if !l.is_finite() || l <= 0.0 { return Err("the solve gave no inductance".into()); }+    Ok((l * 1e9, br.len(), pitch))+}++#[cfg(test)]+/// Validation: a strip (wc cells wide, lc cells long) over a full plane h below, shorted at the far end, port at the near+/// end. Returns (PEEC nH, 2-D field-solve nH from the strip's L per metre over the same length).+pub fn validate_strip(d_mm: f64, wc: usize, lc: usize, h_mm: f64, plane_cells: usize) -> (f64, f64) {+    let w = lc + 2; let hh = plane_cells.max(wc + 2);+    let r0 = (hh - wc) / 2;+    let mut cell = vec![vec![None; w * hh]; 2];+    for r in 0..hh { for c in 0..w { cell[1][r * w + c] = Some(0u16); } }+    for r in r0..r0 + wc { for c in 1..1 + lc { cell[0][r * w + c] = Some(0u16); } }+    let g = Grid { x0: 0.0, y0: 0.0, d: d_mm, w, h: hh, z: vec![h_mm, 0.0], t: 35e-6, cell };+    // the far-end short: one via per strip row at the last strip column+    let vias: Vec<(f64, f64, u16)> = (r0..r0 + wc).map(|r| ((lc as f64 + 0.5) * d_mm, (r as f64 + 0.5) * d_mm, 0u16)).collect();+    let br = g.network(&vias);+    // the port: the near-end strip column against the plane under it, all rows tied by shorts to one node each side+    let mut shorts = Vec::new(); let p = g.node(0, r0, 1); let q = g.node(1, r0, 1);+    for r in r0 + 1..r0 + wc { shorts.push((p, g.node(0, r, 1))); shorts.push((q, g.node(1, r, 1))); }+    let nn = 2 * w * hh;+    let l = solve(nn, &br, &shorts, p, q).unwrap_or(f64::NAN) * 1e9;+    let x = crate::xsec::Xs { w: wc as f64 * d_mm * 1e-3, t: 35e-6, gap: None, below: h_mm * 1e-3, above: None, er: 1.0, mask: None };+    let z = crate::xsec::solve(&x); let lpm = z.z0 / 299_792_458.0;+    (l, lpm * (lc as f64 - 1.0) * d_mm * 1e-3 * 1e9)+}++#[cfg(test)] mod tests {+    use super::*;+    #[test] fn wire_mutual_closed_form() {+        let (l, d) = (0.01, 0.001); let m = m_par(0.0, l, 0.0, l, d);+        let want = 2e-7 * l * ((l / d).asinh() - (1.0 + d * d / (l * l)).sqrt() + d / l);+        assert!((m - want).abs() / want < 1e-12);+    }+    #[test] fn strip_over_plane_matches_the_field_solve() {+        // a 1 mm strip 0.5 mm over a plane, shorted at the far end: within 3 % of the 2-D solve's L per metre × length+        let (pe, fs) = validate_strip(0.5, 2, 20, 0.5, 12); assert!((pe / fs - 1.0).abs() < 0.03, "PEEC {pe} nH, field solve {fs} nH");+    }+}