This is about the inside of the electromagnetic simulator introduced in the previous article. Getting a method-of-moments solver to draw a picture is easy. The hard part is getting to the point where you can say the numbers it produces are right. This article goes through what I checked to get there, and the errors that turned up along the way.
What is being solved
The unknowns are the currents on the conductor surface, and the equation is the mixed-potential integral equation (MPIE), discretised with roof-shaped basis functions spanning two neighbouring cells (rooftop basis functions). At angular frequency ω the matrix has the form
Z(ω) = j(ωL − Ψ/ω) + Zs(ω)·GL is the partial inductance, Ψ the potential coefficient of the charges and Zs the conductor's surface impedance (including the skin effect). L and Ψ do not depend on frequency, so they are assembled once and reused at every frequency. That property is also the key to the speed-ups in the third article.
Qualitative tests miss quantitative errors
The first version passed every qualitative test: "a line passes", "a gap blocks", "a via shorts", "a stub produces a dip". But checking the numbers one by one against independent methods turned up these errors:
| Problem | What was wrong |
|---|---|
| Capacitance up to 40 % too small | The substrate was approximated by the average of air and dielectric. The thinner the substrate, the worse it got, and every resonance frequency was off by about 25 % |
| Numerical integration breaking down | A 6×6-point integration could not capture the sharp peak of the Green's function on a thin substrate (h = 10 µm); the capacitance error reached 83 % |
| A display error | Adding the magnitudes of complex currents made the arrows always point the same way, hiding the reversals of a standing wave |
None of these shows up if you only look at the shape of the S-parameters. Even with every resonance 25 % off, the line still passes and the gap still blocks.
Two things fixed them:
- An exact Green's function: the electrostatic Green's function of a grounded dielectric substrate, evaluated exactly as an infinite series of images. It satisfies both the thin-substrate limit (parallel-plate capacitance) and the interface limit at once
- Closed-form cell integrals: the integrals over cell surfaces are replaced by a closed form written as differences at the corner coordinates, instead of numerical integration. Exact, and the evaluation points dropped from 36 to 4
Only the conductor edges need fine cells
The surface current on a conductor grows as 1/√d towards an edge (d is the distance from the edge). Capturing that sharp variation needs fine cells near the edge.
At first every pattern cell was divided uniformly into M×M. But the singularity only matters very close to the edge. So I changed to keeping the interior coarse and laying layers that get geometrically finer (e, 3e, 7e…) only along conductor edges — the same idea as ADS Momentum's Edge Mesh and Sonnet's edge refinement.
Comparing the characteristic impedance of a one-cell-wide line with the Hammerstad-Jensen closed form:
| Method | Unknowns N | Z₀ error |
|---|---|---|
| No refinement | 18 | 16.5 % |
| Uniform 3×3 | 254 | 5.2 % |
| Uniform 4×4 | 474 | 4.0 % |
| Two layers at the edges only | 186 | 0.2 % |
Fewer unknowns, and more than 20 times more accurate. Since the number of cells is set by the length of the edges rather than the area, larger structures gain more: a 17×17 solid conductor went from 5102 unknowns to 842. The assumption that "more accuracy means finer everywhere" was itself wrong.
What convergence testing showed
That still leaves the question of how to choose the mesh settings: the interior subdivision M, the number of edge layers L and the width e of the innermost layer. I measured the maximum S-parameter error of each setting against a very fine reference mesh (M3, L4, e = 0.02Δ).
| Structure | Setting | N | Max error |
|---|---|---|---|
| Line on the default board, h = 0.8 mm | Old default M1, L2, e = 0.12Δ | 195 | 3.8e-2 |
| Fine edges only M1, L4, e = 0.02Δ | 418 | 7.0e-2 | |
| M2, L4, e = 0.02Δ | 809 | 1.3e-2 | |
| M3, L3, e = 0.03Δ | 1013 | 5.2e-3 | |
| Automatic M2, L3, e = 0.03Δ | 624 | 1.2e-2 | |
| Line on a thin board, h = 0.2 mm | Old default | 195 | 8.2e-2 |
| Automatic M2, L2, e = 0.02Δ | 425 | 2.2e-2 |
Three things came out of it:
- The interior subdivision M matters: the cell size gives 20 cells per wavelength at the design frequency, but only about 11 at the top frequency (6 GHz). With rooftop basis functions that produces a phase error (going from M = 1 to 2 changes things by about 7 %). The earlier advice that "M can normally stay at 1" was wrong
- Refining only the edges is not enough: the setting that kept M1 and refined the edges (N = 418) is actually worse. Convergence in the edge width e is slow, roughly as √e, and tightening only the edges is pointless while the interior stays coarse
- The skin depth is not a mesh scale: it describes the current distribution through the metal's thickness, and that is already in the surface impedance in closed form. What sets the width over which current crowds towards an edge in the plane is the substrate thickness h, the cell size Δ and the metal thickness t
From this I made a rule that sets the mesh automatically (on by default):
- Interior subdivision
M = ⌈20Δ / λg,min⌉(20 cells per wavelength at the top frequency) - Width over which current crowds at an edge
ℓ = min(h, Δ/2), innermost layer widthe = max(ℓ/12, t/2, 0.02Δ)(no finer than half the metal thickness, since the conductor is modelled as a sheet) - The number of layers L is the smallest with
e·(2ᴸ − 1) ≥ ℓ/2(at most 4) - If the unknowns would exceed the limit, coarsen in the order L → M → e and say so
With this rule the error is about a third of the old default's.
The displayed current was not the operating current
To get S-parameters, the solver uses the solution with a voltage on port 1 and port 2 shorted (to extract the Y parameters). The first version drew that solution directly as the surface current.
But a line with a shorted end has a short-circuit resonance every half wavelength. So at certain frequencies the current appeared to jump by 56–79 times. That resonance never appears in the S-parameters, which are taken with 50 Ω terminations. Only the displayed current was showing a state that does not occur in practice.
The anomaly became visible once a feature for comparing currents across frequencies was added. Changing the mesh barely moved the resonance (2.80 → 2.85 GHz), which ruled out the mesh and pointed at the excitation conditions themselves.
Now the two port solutions are combined to draw the current with port 1 driven and port 2 terminated in 50 Ω. On the line, the current magnitude is nearly constant across the band (a ratio of 1.2 between maximum and minimum). Recomputing S11 and S21 from this current alone agrees with the S-parameters to 1e-14.
What is checked
| Item | Result |
|---|---|
| Potential Green's function vs the exact spectral-domain solution | 0.0–0.5 % |
| Effective permittivity vs Hammerstad-Jensen (εr = 2.2–10) | 1–3 % |
| Z₀ of a one-cell-wide line vs Hammerstad-Jensen | 0.2 % |
| Reciprocity (S21 = S12) | 4e-16 |
| Passivity (reflected plus transmitted power no more than 1) | No violations |
| Extreme values (8 cases such as t = 5 mm, σ = 0.01 MS/m, h = 10 µm) | No case breaks down |
| A vertical line vs a horizontal line (rotated 90°) | Agree to 3e-12 dB |
Check, one by one and quantitatively, every quantity that can be compared with a known answer. It is unglamorous, but skip it and you have a tool that merely draws plausible pictures.
The next article is about how this computation was made fast inside the browser. → Solving dense matrices fast in the browser — adaptive frequency sampling and WebGPU