fdm: dielectric-side eps on metal faces + richardson-refine cli paths

sikit/si/CrossSection.h

@@ -45,10 +45,12 @@ struct CrossSection {
     std::vector<DielectricLayer> stack;   // top-to-bottom
     std::vector<Conductor> 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;
 

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]);

sikit/si/TraceImpedance.cpp

@@ -6,6 +6,7 @@
 
 #include <algorithm>
 #include <cmath>
+#include <optional>
 #include <format>
 #include <unordered_map>
 #include <utility>
@@ -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<double> {
+        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;

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);
+}