OCuLink Riser for the Supermicro AIOM Slot

#!/usr/bin/env python3
"""2D finite-difference field solver for edge-coupled differential pairs on the
JLC06161H-3313 stackup (JLCPCB published values, read 2026-10-12 from jlcpcb.com/impedance).

Solves Laplace with layered dielectrics for the odd and even modes, with and without
dielectrics, and reports Zodd, Zeven, Zdiff and Zse. Rectangular conductors (no etch
trapezoid), lossless, quasi-static. Writes design/zsolve.json.

Layer use: L1 sig (ref L2), L2 GND, L3 sig (ref L4), L4 GND, L5 sig (ref L6), L6 GND.
"""
import json, sys
import numpy as np
import scipy.sparse as sp
import scipy.sparse.linalg as spla

C0 = 299792458.0
EPS0 = 8.8541878128e-12
D = 0.0025  # grid step, mm


def solve(regions, conductors, grounds, W, H, excite):
    """regions: list of (x0,x1,y0,y1,er) applied in order (later wins), y up, mm.
    conductors: list of rectangles (x0,x1,y0,y1) for signal 1 and 2.
    grounds: list of rectangles held at 0 V. The outer box is 0 V on top/bottom, Neumann on the sides."""
    nx, ny = int(round(W / D)) + 1, int(round(H / D)) + 1
    # cell permittivity (cells between nodes)
    eps = np.ones((ny - 1, nx - 1))
    xc = (np.arange(nx - 1) + 0.5) * D
    yc = (np.arange(ny - 1) + 0.5) * D
    for (x0, x1, y0, y1, er) in regions:
        mx = (xc >= x0) & (xc <= x1)
        my = (yc >= y0) & (yc <= y1)
        eps[np.ix_(my, mx)] = er
    xn = np.arange(nx) * D
    yn = np.arange(ny) * D
    fixed = np.zeros((ny, nx), bool)
    val = np.zeros((ny, nx))

    def mark(rect, v):
        x0, x1, y0, y1 = rect
        mx = (xn >= x0 - 1e-9) & (xn <= x1 + 1e-9)
        my = (yn >= y0 - 1e-9) & (yn <= y1 + 1e-9)
        fixed[np.ix_(my, mx)] = True
        val[np.ix_(my, mx)] = v

    for g in grounds:
        mark(g, 0.0)
    fixed[0, :] = True; val[0, :] = 0
    fixed[-1, :] = True; val[-1, :] = 0
    for c, v in zip(conductors, excite):
        mark(c, v)

    idx = -np.ones((ny, nx), int)
    free = ~fixed
    idx[free] = np.arange(free.sum())
    n = int(free.sum())
    # padded cell permittivity so out-of-range cells are 0 (Neumann at the sides)
    ep = np.zeros((ny + 1, nx + 1))
    ep[1:-1, 1:-1] = eps
    # edge weights: horizontal edge (j,i)-(j,i+1) uses cells (j-1,i),(j,i) -> ep[j, i+1], ep[j+1, i+1]
    wx = 0.5 * (ep[0:ny, 1:nx] + ep[1:ny + 1, 1:nx])        # shape (ny, nx-1)
    wy = 0.5 * (ep[1:ny, 0:nx] + ep[1:ny, 1:nx + 1])
    diag = np.zeros((ny, nx))
    diag[:, :-1] += wx; diag[:, 1:] += wx
    diag[:-1, :] += wy; diag[1:, :] += wy
    rows, cols, data = [], [], []
    b = np.zeros(n)

    # horizontal couplings
    for (sa, sb, w) in (((slice(None), slice(0, nx - 1)), (slice(None), slice(1, nx)), wx),
                        ((slice(None), slice(1, nx)), (slice(None), slice(0, nx - 1)), wx),
                        ((slice(0, ny - 1), slice(None)), (slice(1, ny), slice(None)), wy),
                        ((slice(1, ny), slice(None)), (slice(0, ny - 1), slice(None)), wy)):
        fa = free[sa]; fb = free[sb]
        ia = idx[sa]; ib = idx[sb]; vb = val[sb]
        m = fa & fb & (w > 0)
        rows.append(ia[m]); cols.append(ib[m]); data.append(-w[m])
        m2 = fa & (~fb) & (w > 0)
        np.add.at(b, ia[m2], w[m2] * vb[m2])
    rows.append(idx[free]); cols.append(idx[free]); data.append(diag[free])
    A = sp.csr_matrix((np.concatenate(data), (np.concatenate(rows), np.concatenate(cols))), shape=(n, n))
    x = spla.spsolve(A, b)
    phi = val.copy()
    phi[free] = x
    # energy per unit length: 0.5 * eps0 * sum eps |E|^2 over cells (cell-centred gradient)
    ex = (phi[:-1, 1:] - phi[:-1, :-1] + phi[1:, 1:] - phi[1:, :-1]) / (2 * D)
    ey = (phi[1:, :-1] - phi[:-1, :-1] + phi[1:, 1:] - phi[:-1, 1:]) / (2 * D)
    Wd = 0.5 * EPS0 * np.sum(eps * (ex ** 2 + ey ** 2)) * (D * D)  # D cancels with 1/D^2: J/m per mm^2/mm^2
    return Wd  # J per metre of line (units consistent: E in V/mm, area in mm^2 -> V^2)


def pair(layer, w, s, er_override=None, with_mask=True):
    """Build geometry for a pair on layer L1, L3 or L5 of JLC06161H-3313."""
    W = 2 * w + s + 1.6
    x1 = W / 2 - s / 2 - w
    c1 = (x1, x1 + w)
    c2 = (W / 2 + s / 2, W / 2 + s / 2 + w)
    if layer == "L1":
        h = 0.0994; t = 0.035; er_pp = 4.1
        top_air = 0.6
        Hh = 0.03 + h + t + 0.03 + top_air
        y_gnd = 0.03  # GND plane L2 slab from 0 to 0.03 (fixed 0)
        y_sig = y_gnd + h
        regions = [(0, W, y_gnd, y_sig, er_pp)]
        if with_mask:
            # conformal solder mask: 1.2 mil over substrate, 0.6 mil over trace (JLC: C1 1.2, C2 0.6, C3 1.2, er 3.8)
            regions.append((0, W, y_sig, y_sig + 0.03048, 3.8))
            for (a, bb) in (c1, c2):
                regions.append((a - 0.01524, bb + 0.01524, y_sig, y_sig + t + 0.01524, 3.8))
        conds = [(c1[0], c1[1], y_sig, y_sig + t), (c2[0], c2[1], y_sig, y_sig + t)]
        grounds = [(0, W, 0, y_gnd)]
        return regions, conds, grounds, W, Hh
    if layer in ("L3", "L5"):
        t = 0.0152
        if layer == "L3":
            h_ref = 0.1088; er_pp = 4.16   # L3 -> L4 through 2116 prepreg
        else:
            h_ref = 0.0994; er_pp = 4.1    # L5 -> L6 through 3313 prepreg
        h_core = 0.55; er_core = 4.6       # the core to the far plane (L2 or L4)
        g = 0.02
        y0 = g
        y_sig = y0 + h_ref - t   # trace bottom sits so its top touches core? copper is on the core side
        # inner copper is etched on the core: trace top face on the core surface, prepreg fills below it
        y_core = y0 + h_ref       # core bottom surface (trace top)
        Hh = y_core + h_core + g
        regions = [(0, W, y0, y_core, er_pp), (0, W, y_core, y_core + h_core, er_core)]
        conds = [(c1[0], c1[1], y_core - t, y_core), (c2[0], c2[1], y_core - t, y_core)]
        grounds = [(0, W, 0, y0), (0, W, Hh - g, Hh)]
        return regions, conds, grounds, W, Hh
    raise ValueError(layer)


def zdiff(layer, w, s):
    reg, conds, gnds, W, H = pair(layer, w, s)
    air = [(r[0], r[1], r[2], r[3], 1.0) for r in reg]
    out = {}
    for mode, ex in (("odd", (1.0, -1.0)), ("even", (1.0, 1.0))):
        Wd = solve(reg, conds, gnds, W, H, ex)
        Wa = solve(air, conds, gnds, W, H, ex)
        # for (+1,-1) or (+1,+1): W = C_mode (per line) * 1V^2  (two lines, 0.5*C*V^2 each)
        Cm = Wd           # J/m: eps0 [F/m] * E^2 [V^2/mm^2] * area [mm^2]
        Ca = Wa
        Z = 1.0 / (C0 * np.sqrt(Cm * Ca))
        out[mode] = {"Z": Z, "eeff": Cm / Ca}
    zd = 2 * out["odd"]["Z"]
    return {"layer": layer, "w": w, "s": s, "Zodd": out["odd"]["Z"], "Zeven": out["even"]["Z"],
            "Zdiff": zd, "eeff_odd": out["odd"]["eeff"], "eeff_even": out["even"]["eeff"]}


if __name__ == "__main__":
    res = []
    cases = [(l, w, s) for l in ("L1", "L3", "L5") for (w, s) in
             [(0.10, 0.15), (0.11, 0.15), (0.12, 0.15), (0.125, 0.15), (0.13, 0.15), (0.14, 0.15),
              (0.10, 0.20), (0.11, 0.20), (0.12, 0.20)]]
    if len(sys.argv) > 1:
        cases = [(sys.argv[1], float(sys.argv[2]), float(sys.argv[3]))]
    for c in cases:
        r = zdiff(*c)
        res.append(r)
        print(f"{r['layer']} w={r['w']:.3f} s={r['s']:.3f}  Zdiff={r['Zdiff']:.1f}  Zodd={r['Zodd']:.1f} Zeven={r['Zeven']:.1f} eeff_odd={r['eeff_odd']:.2f}", flush=True)
    if len(sys.argv) == 1:
        json.dump(res, open("design/zsolve.json", "w"), indent=1)