From 33e21658ea8ca541263e19d9c7fd5659df9f2283 Mon Sep 17 00:00:00 2001 From: UnsignedChad <167274875+UnsignedChad@users.noreply.github.com> Date: Tue, 9 Jun 2026 19:51:28 +0000 Subject: [PATCH] fdm: use dielectric-side eps on metal faces, richardson-refine the cli paths Conductor cells store a placeholder eps_r that was harmonic-meaned into every metal/dielectric face coefficient, wrapping traces in a fake low-eps skin half a cell thick. eps_eff converged ~25% low at the default cell size and the fdm/closed-form engines disagreed 15% on the demo board. Metal faces now take the dielectric cell's eps alone (stencil + Gauss surface stay consistent), production impedance paths use the existing richardson-refined solve, and a regression test pins eps_eff against Hammerstad. Engines now agree to 0.1% on demo_board USB_DP. --- sikit/si/CrossSection.h | 10 +++++---- sikit/si/FdmSolver.cpp | 28 ++++++++++++++++++------ sikit/si/TraceImpedance.cpp | 43 +++++++++++++++++++++++++++---------- sikit/tests/em2d_test.cpp | 42 ++++++++++++++++++++++++++++++++++++ 4 files changed, 101 insertions(+), 22 deletions(-) diff --git a/sikit/si/CrossSection.h b/sikit/si/CrossSection.h index e0b78ea..52ce7a7 100644 --- a/sikit/si/CrossSection.h +++ b/sikit/si/CrossSection.h @@ -45,10 +45,12 @@ struct CrossSection { std::vector stack; // top-to-bottom std::vector conductors; - // Lateral extent of the analysis region. The horizontal walls are - // treated as Neumann (∂V/∂n = 0) boundaries; pick this wide enough - // that field rolls off to ~0 before reaching them (rule of thumb: - // 5–10× the trace width or stack thickness, whichever is larger). + // Lateral extent of the analysis region. All four box walls are + // grounded Dirichlet boundaries (V = 0): solve() iterates interior + // cells only and leaves the boundary ring at its initial 0 V. Pick + // this wide enough that the field rolls off to ~0 before reaching + // them (rule of thumb: 5–10× the trace width or stack thickness, + // whichever is larger), otherwise the grounded box adds capacitance. double y_min = -5e-3; double y_max = 5e-3; diff --git a/sikit/si/FdmSolver.cpp b/sikit/si/FdmSolver.cpp index ee54d6b..d62f9f1 100644 --- a/sikit/si/FdmSolver.cpp +++ b/sikit/si/FdmSolver.cpp @@ -72,7 +72,9 @@ FdmGrid build_grid(const CrossSection& cs, double cell_size_m) { g.is_conductor[k] = 1u; g.conductor_id[k] = cid; g.V[k] = cv; - g.epsilon_r[k] = 1.0; // value doesn't matter inside metal + g.epsilon_r[k] = 1.0; // placeholder; face coefficients never + // read ε from inside metal (see solve() + // and charge_per_length()) } else { g.epsilon_r[k] = cs.epsilon_r_at(yc, zc); } @@ -96,10 +98,20 @@ SolveResult solve(FdmGrid& g, const SolveConfig& cfg) { if (g.is_conductor[k]) continue; const double ec = g.epsilon_r[k]; - const double e_e = face_eps(ec, g.epsilon_r[g.idx(i + 1, j)]); - const double e_w = face_eps(ec, g.epsilon_r[g.idx(i - 1, j)]); - const double e_n = face_eps(ec, g.epsilon_r[g.idx(i, j + 1)]); - const double e_s = face_eps(ec, g.epsilon_r[g.idx(i, j - 1)]); + // Face into metal: the flux path from this cell's centre to + // the metal surface lies entirely in this cell's dielectric, + // so the coefficient is ec alone. Harmonic-meaning against + // the placeholder ε stored inside conductors would wrap every + // trace in a fake low-ε skin half a cell thick (O(h) bias: + // C low, ε_eff low, Z₀ high). + auto fe = [&](std::size_t kn) { + return g.is_conductor[kn] ? ec + : face_eps(ec, g.epsilon_r[kn]); + }; + const double e_e = fe(g.idx(i + 1, j)); + const double e_w = fe(g.idx(i - 1, j)); + const double e_n = fe(g.idx(i, j + 1)); + const double e_s = fe(g.idx(i, j - 1)); const double sum_e = e_e + e_w + e_n + e_s; if (sum_e <= 0.0) continue; @@ -232,8 +244,10 @@ double charge_per_length(const FdmGrid& g, int conductor_id) { auto contribute_face = [&](int ii, int jj) { const std::size_t kk = g.idx(ii, jj); if (g.is_conductor[kk]) return; // facing another metal cell - // Face-averaged ε on the metal/dielectric boundary. - const double e_face = face_eps(g.epsilon_r[k], g.epsilon_r[kk]); + // Dielectric-side ε. Must match the coefficient the solver + // used for this same face (see fe() in solve()): the metal + // cell's stored ε is a placeholder and carries no flux. + const double e_face = g.epsilon_r[kk]; // E field across the face: −(V_neighbour − V_metal)/h, outward. // Flux per unit length = ε · E · (face length) = ε · ΔV. q += kEps0 * e_face * (v_c - g.V[kk]); diff --git a/sikit/si/TraceImpedance.cpp b/sikit/si/TraceImpedance.cpp index e5cda01..8915adc 100644 --- a/sikit/si/TraceImpedance.cpp +++ b/sikit/si/TraceImpedance.cpp @@ -6,6 +6,7 @@ #include #include +#include #include #include #include @@ -240,19 +241,36 @@ double compute_diff_z0_fdm(double trace_width, double spacing, cfg.max_iterations = 100000; const double h = cell_size_for(trace_width, s.copper_thickness); - em2d::FdmGrid g = em2d::build_grid(cs, h); - auto r1 = em2d::solve(g, cfg); - if (!r1.ok) return 0.0; - const double q = em2d::charge_per_length(g, /*trace_p=*/0); - if (q <= 0.0) return 0.0; + // Odd-mode charge on the positive trace at one cell size; nullopt on + // any solver failure. + auto solve_q = [&](const em2d::CrossSection& sec, + double cell) -> std::optional { + em2d::FdmGrid g = em2d::build_grid(sec, cell); + if (!em2d::solve(g, cfg).ok) return std::nullopt; + const double q = em2d::charge_per_length(g, /*trace_p=*/0); + if (q <= 0.0) return std::nullopt; + return q; + }; auto cs_air = cs; for (auto& d : cs_air.stack) d.epsilon_r = 1.0; - em2d::FdmGrid g_air = em2d::build_grid(cs_air, h); - auto r2 = em2d::solve(g_air, cfg); - if (!r2.ok) return 0.0; - const double q_air = em2d::charge_per_length(g_air, /*trace_p=*/0); - if (q_air <= 0.0) return 0.0; + + const auto qf = solve_q(cs, h); + const auto qaf = solve_q(cs_air, h); + if (!qf || !qaf) return 0.0; + + // First-order Richardson against a 2h solve, same scheme as + // compute_z0_refined. Falls back to the fine-only values if either + // coarse solve fails or the extrapolation goes non-physical. + double q = *qf; + double q_air = *qaf; + const auto qc = solve_q(cs, 2.0 * h); + const auto qac = solve_q(cs_air, 2.0 * h); + if (qc && qac) { + const double qe = 2.0 * *qf - *qc; + const double qae = 2.0 * *qaf - *qac; + if (qe > 0.0 && qae > 0.0) { q = qe; q_air = qae; } + } // Excitation is V_p = +0.5, V_n = −0.5 (V_diff = 1). For pure odd mode // the per-conductor odd-mode capacitance is C_odd = Q / (V_p) = 2·Q. @@ -275,7 +293,10 @@ SegmentImpedance compute_one_fdm(double trace_width, int layer_ordinal, cfg.tolerance = 5e-6; cfg.max_iterations = 100000; const double h = cell_size_for(trace_width, s.copper_thickness); - auto out = em2d::compute_z0(cs, 0, 1, h, cfg); + // Refined = solve at h and 2h, Richardson-extrapolate C and C_air. + // Cancels the leading O(h) discretisation error (cell-rect edge + // fattening, metal-surface offset) that a single solve carries. + auto out = em2d::compute_z0_refined(cs, 0, 1, h, cfg); if (!out.ok) { r.in_valid_range = false; return r; diff --git a/sikit/tests/em2d_test.cpp b/sikit/tests/em2d_test.cpp index 8901b2a..38e1693 100644 --- a/sikit/tests/em2d_test.cpp +++ b/sikit/tests/em2d_test.cpp @@ -124,3 +124,45 @@ TEST_CASE("em2d: microstrip Z₀ within 15% of Wadell closed-form", "[em2d]") { REQUIRE(r.eps_eff > 2.0); // somewhere between air (1) and FR-4 (4.4) REQUIRE(r.eps_eff < 4.4); } + +TEST_CASE("em2d: microstrip eps_eff matches Hammerstad closed form", "[em2d]") { + // Regression guard for the metal-face ε bug: harmonic-meaning the + // placeholder ε stored inside conductor cells into surface face + // coefficients wrapped every trace in a fake low-ε skin half a cell + // thick, dragging eps_eff for this geometry down to ~2.4 and Z₀ up + // ~15%. Hammerstad for w/h = 1.84, εr = 4.4 gives eps_eff ≈ 3.32; + // the solver must land on the physical side of 3.0. + CrossSection cs; + cs.stack.push_back({1.524e-3, 4.4, 0, "FR4"}); + + Conductor trace; + trace.id = 0; + trace.y_center = 0; + trace.z_top = -35e-6; // surface trace: sits in air, bottom on FR-4 + trace.width = 2.8e-3; + trace.thickness = 35e-6; + trace.voltage = 1.0; + + Conductor gnd; + gnd.id = 1; + gnd.y_center = 0; + gnd.z_top = 1.524e-3; + gnd.width = 20e-3; + gnd.thickness = 35e-6; + gnd.voltage = 0.0; + + cs.conductors.push_back(trace); + cs.conductors.push_back(gnd); + cs.y_min = -10e-3; + cs.y_max = 10e-3; + cs.air_above = 5e-3; + + SolveConfig cfg; + cfg.tolerance = 5e-6; + cfg.max_iterations = 80000; + auto r = compute_z0(cs, 0, 1, 100e-6, cfg); + + REQUIRE(r.ok); + REQUIRE(r.eps_eff > 3.0); + REQUIRE(r.eps_eff < 3.7); +}