emikit: openEMS cross-check scripts (WIP)

emikit/openems/README.md

@@ -0,0 +1,34 @@
+# openEMS cross-check scripts
+
+Full-wave FDTD against emikit's closed form on the same single-trace
+fixture. If the two agree to a few dB on the geometry emikit models,
+the chamber-vs-emikit gap on the TI EVM is scope (multi-net, activity
+factor, cable CM), not math.
+
+## status
+
+scripts run, openEMS builds and produces port + near-field data, but
+the broadband Gaussian + -40 dB auto-terminate stops at ~2 ns -- too
+short for a clean 30 MHz spectrum. needs more iteration before the
+numbers are usable.
+
+files:
+
+* `openems_microstrip.py` -- 30 mm trace, 0.15 mm wide, 0.2 mm above
+  solid GND, lumped ports + NF2FF at 10 m
+* `compare_openems_emikit.py` -- runs emikit's `loop_e_field` over the
+  FDTD-measured I_port spectrum
+
+run on a workstation with openEMS Python bindings:
+
+    source ~/opt/openEMS/venv/bin/activate
+    python3 openems_microstrip.py
+    python3 compare_openems_emikit.py
+
+## todo
+
+- bump `NrTS`, disable `EndCriteria`, or run CW per frequency instead
+  of broadband Gaussian
+- try `AddMSLPort` instead of `AddLumpedPort` (less sensitivity near
+  the port terminals)
+- mesh capped under 2M cells; tighter spacing blows past 30 GB

emikit/openems/compare_openems_emikit.py

@@ -0,0 +1,57 @@
+#!/usr/bin/env python3
+"""feeds openEMS's measured port-current spectrum into emikit's
+loop_e_field and compares to the full-wave far-field."""
+import math
+import h5py
+import numpy as np
+
+ETA0 = 376.730313668
+C0   = 2.99792458e8
+
+# Geometry (must match openems_run.py).
+TRACE_LEN_M    = 30.0e-3
+LOOP_HEIGHT_M  = 0.2e-3
+OBS_DIST_M     = 10.0
+LOOP_AREA      = TRACE_LEN_M * LOOP_HEIGHT_M
+
+
+def emikit_loop_e_field(i_a, freq_hz, distance_m, area_m2):
+    """Ott Eq 11-2 -- same as emikit's LoopEmissions::loop_e_field."""
+    return (ETA0 * math.pi * i_a * area_m2 * freq_hz ** 2) / (C0 ** 2 * distance_m)
+
+
+def main():
+    with h5py.File("/home/chad/openems-runs/emikit_microstrip_air_result.h5",
+                   "r") as h:
+        freqs   = h["freq_hz"][...]
+        i_port  = h["i_port_a"][...]
+        e_full  = h["e_v_m"][...]
+        e_full_dbuv = h["e_dbuv_m"][...]
+
+    print(f"{'freq_MHz':>9} {'|I|mA':>8} {'E_FW_dBuV':>11} "
+          f"{'E_emikit_dBuV':>14} {'delta_dB':>10}")
+    print("-" * 60)
+
+    e_emikit_all = []
+    deltas = []
+    for f, ip, efw_v, efw_db in zip(freqs, i_port, e_full, e_full_dbuv):
+        e_emikit_v = emikit_loop_e_field(abs(ip), f, OBS_DIST_M, LOOP_AREA)
+        if e_emikit_v <= 0:
+            e_emikit_db = -1000.0
+        else:
+            e_emikit_db = 20.0 * math.log10(e_emikit_v * 1e6)
+        e_emikit_all.append(e_emikit_db)
+        delta = e_emikit_db - efw_db
+        deltas.append(delta)
+        print(f"{f / 1e6:9.1f} {abs(ip) * 1e3:8.4f} {efw_db:11.2f} "
+              f"{e_emikit_db:14.2f} {delta:+10.2f}")
+
+    valid = [d for d in deltas if abs(d) < 200]
+    if valid:
+        print(f"\nmean |delta| = {np.mean(np.abs(valid)):.2f} dB")
+        print(f"max  |delta| = {np.max (np.abs(valid)):.2f} dB")
+        print(f"rms       = {np.sqrt(np.mean(np.array(valid) ** 2)):.2f} dB")
+
+
+if __name__ == "__main__":
+    main()

emikit/openems/openems_microstrip.py

@@ -0,0 +1,104 @@
+#!/usr/bin/env python3
+"""Full-fixture FDTD with sane mesh grading."""
+import os
+import numpy as np
+import h5py
+from openEMS import openEMS
+from CSXCAD import ContinuousStructure
+
+sim = os.path.expanduser("~/openems-runs/full")
+os.makedirs(sim, exist_ok=True)
+
+TRACE_LEN = 30.0
+TRACE_W   = 0.15
+TRACE_T   = 0.035
+SUB_H     = 0.2
+GND_HALF  = 30.0
+AIR       = 80.0
+
+FDTD = openEMS(EndCriteria=1e-4, NrTS=60000)
+FDTD.SetGaussExcite(500e6, 500e6)
+FDTD.SetBoundaryCond(['MUR'] * 6)
+CSX = ContinuousStructure()
+FDTD.SetCSX(CSX)
+CSX.GetGrid().SetDeltaUnit(1e-3)
+
+gnd = CSX.AddMetal("gnd")
+gnd.AddBox([-GND_HALF, -GND_HALF, 0.0], [GND_HALF, GND_HALF, 0.0])
+
+trace = CSX.AddMetal("trace")
+trace.AddBox([-TRACE_LEN / 2, -TRACE_W / 2, SUB_H],
+             [ TRACE_LEN / 2,  TRACE_W / 2, SUB_H + TRACE_T])
+
+FDTD.AddLumpedPort(1, 50.0,
+                   [-TRACE_LEN / 2, -TRACE_W / 2, 0.0],
+                   [-TRACE_LEN / 2,  TRACE_W / 2, SUB_H],
+                   'z', excite=True)
+FDTD.AddLumpedPort(2, 50.0,
+                   [ TRACE_LEN / 2, -TRACE_W / 2, 0.0],
+                   [ TRACE_LEN / 2,  TRACE_W / 2, SUB_H],
+                   'z', excite=False)
+
+# Explicit mesh lines near the trace + coarse at boundary, then smooth
+# with a generous max_res (2 mm) so the mesh GROWS away from the fine
+# lines instead of being held at 0.05 mm everywhere.
+mesh = CSX.GetGrid()
+
+# X: coarse end-of-domain + fine ends-of-trace + middle filler.
+x_lines = [-GND_HALF - AIR, -GND_HALF, -TRACE_LEN/2 - 1.0,
+           -TRACE_LEN/2 - 0.2, -TRACE_LEN/2,
+           -TRACE_LEN/4, 0.0, TRACE_LEN/4,
+           TRACE_LEN/2, TRACE_LEN/2 + 0.2,
+           TRACE_LEN/2 + 1.0, GND_HALF, GND_HALF + AIR]
+mesh.AddLine('x', x_lines)
+mesh.SmoothMeshLines('x', 2.0, ratio=1.4)
+
+# Y: tight cluster around trace, then coarse out.
+y_lines = [-GND_HALF - AIR, -GND_HALF, -2.0, -0.5, -TRACE_W,
+           -TRACE_W/2, -TRACE_W/4, 0.0, TRACE_W/4, TRACE_W/2,
+           TRACE_W, 0.5, 2.0, GND_HALF, GND_HALF + AIR]
+mesh.AddLine('y', y_lines)
+mesh.SmoothMeshLines('y', 2.0, ratio=1.4)
+
+# Z: fine through dielectric+trace, coarse out.
+z_lines = [-AIR, -10.0, -1.0, 0.0, SUB_H/4, SUB_H/2, 3*SUB_H/4,
+           SUB_H, SUB_H + TRACE_T/2, SUB_H + TRACE_T,
+           SUB_H + 0.2, 1.0, 10.0, SUB_H + AIR]
+mesh.AddLine('z', z_lines)
+mesh.SmoothMeshLines('z', 2.0, ratio=1.4)
+
+nx, ny, nz = (len(mesh.GetLines(a)) for a in 'xyz')
+cells = nx * ny * nz
+print(f"Mesh: {nx} x {ny} x {nz} = {cells:,} cells", flush=True)
+assert cells < 30_000_000, f"mesh too big ({cells:,}), tune ratios"
+
+nf2ff = FDTD.CreateNF2FFBox()
+
+FDTD.Run(sim, verbose=2, cleanup=False)
+
+freqs = np.linspace(30e6, 1.0e9, 40)
+ff = nf2ff.CalcNF2FF(sim, freqs, np.array([90.0]), np.array([0.0]),
+                     radius=10.0)
+e_ff = np.abs(ff.E_norm[0])[:, 0, 0]
+
+ut = np.loadtxt(os.path.join(sim, 'port_ut_1'))
+it = np.loadtxt(os.path.join(sim, 'port_it_1'))
+dt = ut[1, 0] - ut[0, 0]
+t = ut[:, 0]
+i_t = it[:, 1]
+i_f = np.array([np.sum(i_t * np.exp(-2j * np.pi * f * t)) * dt
+                 for f in freqs])
+i_port = np.abs(i_f)
+
+out = os.path.expanduser("~/openems-runs/full_result.h5")
+with h5py.File(out, "w") as h:
+    h.create_dataset("freq_hz",  data=freqs)
+    h.create_dataset("i_port_a", data=i_port)
+    h.create_dataset("e_v_m",    data=e_ff)
+    h.create_dataset("e_dbuv_m", data=20.0 * np.log10(e_ff * 1e6))
+
+print(f"\n{'freq_MHz':>9} {'|I|mA':>9} {'E_dBuV/m':>10}", flush=True)
+for k in range(len(freqs)):
+    print(f"{freqs[k] / 1e6:9.1f} {i_port[k] * 1e3:9.4f} "
+          f"{20.0 * np.log10(e_ff[k] * 1e6):10.2f}", flush=True)
+print(f"\nwrote {out}", flush=True)