diff --git a/emikit/openems/README.md b/emikit/openems/README.md new file mode 100644 index 0000000..de1fa6c --- /dev/null +++ b/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 diff --git a/emikit/openems/compare_openems_emikit.py b/emikit/openems/compare_openems_emikit.py new file mode 100755 index 0000000..5d6b9ea --- /dev/null +++ b/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() diff --git a/emikit/openems/openems_microstrip.py b/emikit/openems/openems_microstrip.py new file mode 100755 index 0000000..8458117 --- /dev/null +++ b/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)