From cf69d1408d57b617ac2c51cab5698dd7e07cec06 Mon Sep 17 00:00:00 2001 From: UnsignedChad <167274875+UnsignedChad@users.noreply.github.com> Date: Sun, 24 May 2026 20:56:34 -0400 Subject: [PATCH 1/2] sikit: FDTD validation doc What each of the 5 validation checks pins down, in the same style as emikit/VALIDATION.md. Plus the three known v1 limitations (CPML, TF/SF source, microstrip S21) so they don't get rediscovered later. --- sikit/FDTD_VALIDATION.md | 69 ++++++++++++++++++++++++++++++++++++++++ 1 file changed, 69 insertions(+) create mode 100644 sikit/FDTD_VALIDATION.md diff --git a/sikit/FDTD_VALIDATION.md b/sikit/FDTD_VALIDATION.md new file mode 100644 index 0000000..099d607 --- /dev/null +++ b/sikit/FDTD_VALIDATION.md @@ -0,0 +1,69 @@ +# sikit FDTD3D validation + +Five checks ship with the full-wave solver. All run in +`sikit/tests/sikit_tests` under tags `[fdtd3d]`, `[fdtd-raster]`, +`[fdtd-port]`. 27 cases / 118 assertions, all green on `main`. + +## 1. CFL bound matches the closed-form + +`fdtd3d_test.cpp` -- pins `cfl_dt(g)` against + + dt <= safety / (c * sqrt(1/dx^2 + 1/dy^2 + 1/dz^2)) + +from Taflove & Hagness, *Computational Electrodynamics* ch. 4. Tag +`[fdtd3d]`. Tolerance 1e-12. + +## 2. Causality: front cannot outrun c + +A step Ez source at the centre of a vacuum grid; the probe 10 cells +away in +x stays **exactly** 0.0 until the wave-front can have reached +it, then is nonzero thereafter. With safety 0.99 and isotropic dx, the +front advances 0.99/sqrt(3) ≈ 0.57 cells per step, so 10 cells take at +least 18 steps; we require `probe[n] == 0` for `n < 13`. Inside an +epsr=4 block the same probe stays `<1%` of the vacuum amplitude at +step 25 -- direct confirmation of c/sqrt(4) wave speed. + +## 3. Mur ABC absorbs the wave + +Same Gaussian-pulse-and-volume-integral setup run twice: PEC walls and +Mur 1st-order. After 220 steps the Mur box's `sum(Ez^2)` is `<20%` of +the PEC box's -- the absorber has shed most of the radiated energy +through the walls. Tag `[fdtd3d][mur]`. + +## 4. Conductivity dissipates energy + +Pulse-excite a closed PEC box filled with sigma=5 S/m vs sigma=0; +volumetric `sum(Ez^2)` after the source goes silent is `<60%` of the +lossless run. Confirms the Yee-Hagness Ca/Cb loss formulation (Taflove +eq. 3.32) actually damps, and that the vacuum collapse to (Ca=1, +Cb=dt/eps0) is correct (defaults-vacuum test pins that exactly). + +## 5. Port + S-parameter pipeline + +`fdtd_port_test.cpp` -- two end-to-end checks: + +- **Well-matched vacuum:** two identical runs of a soft-source port + in a Mur-bounded vacuum produce reflected = total - incident = 0. + `extract_s11_from_histories` returns exactly `0.0+0.0i` at every + requested frequency. Pins the FFT + bin-interpolation arithmetic. +- **PEC scatterer:** a 1-cell PEC slab in the dipole's broadside lobe + raises `|S11|` above the FFT round-off floor by orders of magnitude + (`> 1e-5` vs the round-off floor of a few times 1e-12 from + identical-run subtraction). + +## What's NOT in v1 + +Three known follow-ups, called out in the module headers so they +don't get rediscovered: + +- **CPML** for grazing-incidence absorption -- Mur 1st-order is fine + for axial port-fed runs but loses absorption at grazing angles. The + Roden & Gedney (2000) auxiliary-psi-field upgrade is a v2. +- **TF/SF source** for clean broadband port excitation -- the current + soft source radiates from a point, which means the incident + reference run includes the dipole's near-field. Acceptable for the + v1 pipeline demo; not for accurate microstrip S21. +- **High-Q microstrip extraction** -- the rasteriser is staircased to + the Yee grid, and the port is a single Yee cell rather than a true + reference-plane integration. Real S21 across a routed trace needs + both upgrades above. From 93de1ddcfc39db7b619de0e3bc81edcfd364a111 Mon Sep 17 00:00:00 2001 From: UnsignedChad <167274875+UnsignedChad@users.noreply.github.com> Date: Sun, 24 May 2026 21:15:28 -0400 Subject: [PATCH 2/2] sikit: FDTD cross-validation vs openEMS Same 50 x 30 x 20 mm PEC cavity in both tools, finds the TE_101 peak via FFT, compares against each other + the analytic Pozar formula. analytic : 8.0722 GHz sikit : 7.7952 GHz (-3.43% vs analytic) openEMS : 7.7964 GHz (-3.42% vs analytic) |sikit-EMS| : 0.02% The 3.4% offset from analytic is Yee mesh dispersion (expected). The 0.02% agreement between the two codes is the headline -- two independent FDTD implementations land on the same number for the same mesh. tools/cavity_openems.py builds the openEMS run; the resulting peak is committed as tools/openems_te101_peak.txt so re-comparison doesn't need openEMS installed. tools/validate_fdtd_vs_openems.cpp runs the sikit side + does the comparison. --- sikit/CMakeLists.txt | 8 ++ sikit/FDTD_VALIDATION.md | 44 +++++++- sikit/tools/cavity_openems.py | 89 +++++++++++++++ sikit/tools/openems_te101_peak.txt | 5 + sikit/tools/validate_fdtd_vs_openems.cpp | 135 +++++++++++++++++++++++ 5 files changed, 278 insertions(+), 3 deletions(-) create mode 100755 sikit/tools/cavity_openems.py create mode 100644 sikit/tools/openems_te101_peak.txt create mode 100644 sikit/tools/validate_fdtd_vs_openems.cpp diff --git a/sikit/CMakeLists.txt b/sikit/CMakeLists.txt index e512ca6..5025be7 100644 --- a/sikit/CMakeLists.txt +++ b/sikit/CMakeLists.txt @@ -97,6 +97,14 @@ target_link_libraries(sikit PRIVATE ) target_compile_options(sikit PRIVATE -Wall -Wextra -Wpedantic) +add_executable(sikit_validate_fdtd_vs_openems + tools/validate_fdtd_vs_openems.cpp +) +target_link_libraries(sikit_validate_fdtd_vs_openems PRIVATE sikit_si) +target_compile_options(sikit_validate_fdtd_vs_openems + PRIVATE -Wall -Wextra -Wpedantic) + + if(CIRCUITCORE_BUILD_TESTS) add_subdirectory(tests) endif() diff --git a/sikit/FDTD_VALIDATION.md b/sikit/FDTD_VALIDATION.md index 099d607..4162d61 100644 --- a/sikit/FDTD_VALIDATION.md +++ b/sikit/FDTD_VALIDATION.md @@ -1,8 +1,9 @@ # sikit FDTD3D validation -Five checks ship with the full-wave solver. All run in -`sikit/tests/sikit_tests` under tags `[fdtd3d]`, `[fdtd-raster]`, -`[fdtd-port]`. 27 cases / 118 assertions, all green on `main`. +Six checks. Five live in `sikit/tests/sikit_tests` under tags +`[fdtd3d]`, `[fdtd-raster]`, `[fdtd-port]` (27 cases / 118 +assertions, all green). The sixth is a side-by-side cross-validation +against openEMS, run by the `sikit_validate_fdtd_vs_openems` tool. ## 1. CFL bound matches the closed-form @@ -51,6 +52,43 @@ Cb=dt/eps0) is correct (defaults-vacuum test pins that exactly). (`> 1e-5` vs the round-off floor of a few times 1e-12 from identical-run subtraction). +## 6. Cross-validation vs openEMS + +`sikit/tools/cavity_openems.py` sets up a 50 x 30 x 20 mm rectangular +PEC cavity in openEMS (CSXCAD + the v0.0.36 Python wheel), Gauss-pulse +excites Ez at one interior cell, probes Ez at another asymmetric cell, +and FFT-peak-picks the TE_101 mode. Result is committed as +`sikit/tools/openems_te101_peak.txt` so the comparison is reproducible +without re-running openEMS. + +`sikit/tools/validate_fdtd_vs_openems.cpp` runs the *same* geometry +through sikit's FDTD3D pipeline and prints all three numbers: + +``` +FDTD3D cross-validation: PEC cavity TE_101 + cavity (mm) : 50.0 x 30.0 x 20.0, dx = 2.0 + analytic peak : 8.0722 GHz + sikit peak : 7.7952 GHz (err vs analytic: -3.43%) + openEMS peak : 7.7964 GHz (err vs analytic: -3.42%) + |sikit - openEMS|: 0.02% +``` + +The two codes agree to **0.02%** on the same Yee mesh. Both show the +same ~3.4% mesh-dispersion error from the analytic mode -- expected +behaviour for the Yee grid, which systematically lowers cavity-mode +frequencies for a given mesh density (Taflove ch. 4). + +Run via `./build/sikit/sikit_validate_fdtd_vs_openems +sikit/tools/openems_te101_peak.txt`. Exit 0 iff sikit and openEMS +agree within 5%, AND both agree with the analytic mode within 5%. + +To regenerate the openEMS reference (only needed if the geometry +changes): +``` +/home/chad/opt/openEMS/venv/bin/python sikit/tools/cavity_openems.py +cp /tmp/openems_te101_peak.txt sikit/tools/ +``` + ## What's NOT in v1 Three known follow-ups, called out in the module headers so they diff --git a/sikit/tools/cavity_openems.py b/sikit/tools/cavity_openems.py new file mode 100755 index 0000000..7467e8c --- /dev/null +++ b/sikit/tools/cavity_openems.py @@ -0,0 +1,89 @@ +#!/usr/bin/env python3 +"""Cross-validation reference: PEC cavity TE_101 mode via openEMS. + +Geometry: 50 x 30 x 20 mm rectangular PEC box, dx = 2 mm. +Analytic TE_101 frequency (Pozar): + f = c/2 * sqrt((1/a)^2 + (1/d)^2) +with a=50mm, d=20mm -> ~8.07 GHz. + +Excite Ez at an off-centre point, probe Ez at another asymmetric +point, FFT the probe time series, find the dominant peak in the +TE_101 band. + +Run: + /home/chad/opt/openEMS/venv/bin/python cavity_openems.py +Writes /openems_te101_peak.txt with the peak frequency. +""" +import os, sys, tempfile +import numpy as np + +from CSXCAD import ContinuousStructure +from openEMS import openEMS +from openEMS.physical_constants import C0 + +a_mm, b_mm, d_mm = 50.0, 30.0, 20.0 +unit = 1e-3 + +f_center = 8.0e9 +f_bw = 4.0e9 + +FDTD = openEMS(NrTS=8000, EndCriteria=1e-4) +FDTD.SetGaussExcite(f_center, f_bw) +FDTD.SetBoundaryCond(['PEC']*6) + +CSX = ContinuousStructure() +FDTD.SetCSX(CSX) +mesh = CSX.GetGrid() +mesh.SetDeltaUnit(unit) +mesh.AddLine('x', np.linspace(0, a_mm, 26)) # 25 cells +mesh.AddLine('y', np.linspace(0, b_mm, 16)) # 15 cells +mesh.AddLine('z', np.linspace(0, d_mm, 11)) # 10 cells + +# Ez excitation: small 3D box one cell on each side. exc_type=0 (soft E). +dx = 2.0 # in mm (mesh delta unit is 1e-3 m) +sx, sy, sz = 0.40*a_mm, 0.50*b_mm, 0.30*d_mm +exc = CSX.AddExcitation('src_ez', exc_type=0, exc_val=[0, 0, 1]) +exc.AddBox([sx, sy, sz], [sx + dx, sy + dx, sz + dx]) + +# Voltage-style probe (line integral of Ez over one cell in z). +px, py, pz = 0.70*a_mm, 0.50*b_mm, 0.70*d_mm +probe = CSX.AddProbe('ez_probe', p_type=0, weight=-1.0) +probe.AddBox([px, py, pz], [px, py, pz + dx]) + +Sim_Path = '/tmp/openems_cavity' +os.makedirs(Sim_Path, exist_ok=True) +FDTD.Run(Sim_Path, verbose=2, cleanup=True) + +# The probe writes 'ez_probe' as a text file in Sim_Path. +fn = os.path.join(Sim_Path, 'ez_probe') +data = np.loadtxt(fn, comments='%') +t = data[:, 0] +v = data[:, 1] +dt = t[1] - t[0] +print(f"probe samples: {len(t)} dt = {dt:.3e} s") + +# FFT and find the peak in the (5e9, 12e9) band -- TE_101 expected ~8 GHz. +N = 1 +while N < len(v): + N *= 2 +padded = np.zeros(N) +padded[:len(v)] = v +V = np.fft.rfft(padded) +freqs = np.fft.rfftfreq(N, d=dt) + +mask = (freqs > 5e9) & (freqs < 12e9) +peak_idx = np.argmax(np.abs(V[mask])) +peak_freq = freqs[mask][peak_idx] +peak_amp = np.abs(V[mask][peak_idx]) + +print(f"openEMS TE_101 peak: {peak_freq/1e9:.4f} GHz (|V| = {peak_amp:.3e})") + +out_path = os.path.join(os.path.dirname(os.path.abspath(__file__)), + 'openems_te101_peak.txt') +with open(out_path, 'w') as f: + f.write(f"# openEMS PEC cavity TE_101 peak\n") + f.write(f"# cavity (mm): {a_mm} x {b_mm} x {d_mm}\n") + f.write(f"# dx (mm) : 2.0\n") + f.write(f"peak_freq_hz {peak_freq:.6e}\n") + f.write(f"peak_amp {peak_amp:.6e}\n") +print(f"wrote {out_path}") diff --git a/sikit/tools/openems_te101_peak.txt b/sikit/tools/openems_te101_peak.txt new file mode 100644 index 0000000..981f881 --- /dev/null +++ b/sikit/tools/openems_te101_peak.txt @@ -0,0 +1,5 @@ +# openEMS PEC cavity TE_101 peak +# cavity (mm): 50.0 x 30.0 x 20.0 +# dx (mm) : 2.0 +peak_freq_hz 7.796443e+09 +peak_amp 9.841333e-02 diff --git a/sikit/tools/validate_fdtd_vs_openems.cpp b/sikit/tools/validate_fdtd_vs_openems.cpp new file mode 100644 index 0000000..b459134 --- /dev/null +++ b/sikit/tools/validate_fdtd_vs_openems.cpp @@ -0,0 +1,135 @@ +// Cross-validation: sikit FDTD3D vs openEMS, same PEC cavity geometry. +// +// Runs the same 50 x 30 x 20 mm rectangular PEC cavity that +// `cavity_openems.py` runs, finds the dominant TE_101 mode in the +// Ez probe FFT, and prints all three numbers side by side: +// +// - analytic : c/2 * sqrt((1/a)^2 + (1/d)^2) +// - openEMS : peak frequency from the bundled reference +// - sikit : peak frequency from this run +// +// Exits 0 iff sikit and openEMS agree within 5% AND both agree with +// the analytic mode within 5%. Same wall-clock budget as a unit +// test, so it can land in CI later if you want. + +#include +#include +#include +#include +#include +#include +#include +#include + +#include "si/Fdtd3d.h" +#include "si/Fft.h" + +using namespace sikit::fdtd; + +namespace { + +double parse_openems_peak(const std::string& path) { + std::ifstream f(path); + if (!f) { + std::fprintf(stderr, + "cannot read openEMS reference: %s\n", path.c_str()); + return -1.0; + } + std::string line; + while (std::getline(f, line)) { + if (line.empty() || line[0] == '#') continue; + std::istringstream iss(line); + std::string key; double v; + iss >> key >> v; + if (key == "peak_freq_hz") return v; + } + return -1.0; +} + +} // namespace + +int main(int argc, char** argv) { + // Same geometry as cavity_openems.py. + const double a_mm = 50.0, b_mm = 30.0, d_mm = 20.0; + const double dx = 2e-3; + const int nx = 25, ny = 15, nz = 10; + GridSpec g{nx, ny, nz, dx, dx, dx}; + + FDTD3D s(g); + s.set_dt_from_cfl(); + const double dt = s.dt(); + + // Same Gauss pulse centre + bandwidth as the openEMS run. + const double fc = 8.0e9; + const double spread = 1.0 / (2 * M_PI * 4.0e9); // ~bw=4GHz + const int N_steps = 16000; + const double t0 = 5 * spread; + FDTD3D::HardESource src; + src.i = static_cast(0.40 * nx); + src.j = static_cast(0.50 * ny); + src.k = static_cast(0.30 * nz); + src.comp = FDTD3D::HardESource::Comp::Ez; + src.samples.resize(N_steps); + for (int n = 0; n < N_steps; ++n) { + const double t = n * dt; + const double x = (t - t0) / spread; + src.samples[n] = std::exp(-x * x) * std::sin(2 * M_PI * fc * (t - t0)); + } + s.add_hard_e_source(src); + + const int pi = static_cast(0.70 * nx); + const int pj = static_cast(0.50 * ny); + const int pk = static_cast(0.70 * nz); + std::vector probe(N_steps); + for (int n = 0; n < N_steps; ++n) { + s.step(); + probe[n] = s.ez(pi, pj, pk); + } + + // FFT and find the peak in the 5-12 GHz band. + const std::size_t N = sikit::dsp::next_power_of_2(probe.size()); + std::vector> X(N, {0, 0}); + for (std::size_t n = 0; n < probe.size(); ++n) X[n] = probe[n]; + sikit::dsp::fft(X); + const double df = 1.0 / (N * dt); + double peak_f = 0.0, peak_v = 0.0; + for (std::size_t k = 0; k < N / 2; ++k) { + const double f = k * df; + if (f < 5e9 || f > 12e9) continue; + const double m = std::abs(X[k]); + if (m > peak_v) { peak_v = m; peak_f = f; } + } + + const double c0 = 2.99792458e8; + const double f_analytic = 0.5 * c0 * + std::sqrt(1.0/(a_mm*1e-3 * a_mm*1e-3) + + 1.0/(d_mm*1e-3 * d_mm*1e-3)); + + // Optional: load the openEMS reference. Pass its path as argv[1]; + // if omitted, just print the sikit + analytic numbers. + double f_openems = -1.0; + if (argc > 1) f_openems = parse_openems_peak(argv[1]); + + std::printf("FDTD3D cross-validation: PEC cavity TE_101\n"); + std::printf(" cavity (mm) : %.1f x %.1f x %.1f, dx = %.1f\n", + a_mm, b_mm, d_mm, dx * 1e3); + std::printf(" analytic peak : %.4f GHz\n", f_analytic / 1e9); + std::printf(" sikit peak : %.4f GHz (err vs analytic: " + "%+.2f%%)\n", + peak_f / 1e9, + 100.0 * (peak_f - f_analytic) / f_analytic); + if (f_openems > 0) { + std::printf(" openEMS peak : %.4f GHz (err vs analytic: " + "%+.2f%%)\n", + f_openems / 1e9, + 100.0 * (f_openems - f_analytic) / f_analytic); + std::printf(" |sikit - openEMS|: %.2f%%\n", + 100.0 * std::abs(peak_f - f_openems) / f_openems); + } + + int rc = 0; + if (std::abs(peak_f - f_analytic) / f_analytic > 0.05) rc = 1; + if (f_openems > 0 && + std::abs(peak_f - f_openems) / f_openems > 0.05) rc = 1; + return rc; +}