.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "howto/plot_mesh_convergence.py" .. LINE NUMBERS ARE GIVEN BELOW. .. only:: html .. note:: :class: sphx-glr-download-link-note :ref:`Go to the end ` to download the full example code. .. rst-class:: sphx-glr-example-title .. _sphx_glr_howto_plot_mesh_convergence.py: Mesh convergence: a ladder around your simulation ================================================= Nothing inside a single run tells you how far its result is from the converged one. The instrument that does is a *ladder*: the same model on successively finer meshes, one number read off each rung, and a stop rule for the number. Commercial suites hide this loop behind an "adaptive mesh refinement" button; this page is the loop, written out, so it wraps around whatever you are simulating. Three blocks. **Before** your simulation: one function that turns a rung number into a ``MeshControl``. **Your simulation**: unchanged, except that it takes that ``MeshControl`` instead of a hard-coded one. **After** it: record the quantity you design with, compare it with the previous rung, stop when the change is below your tolerance twice in a row. The first half of the page does this for S-parameters, the second for a resonance frequency; both blocks copy into any script. .. GENERATED FROM PYTHON SOURCE LINES 21-33 .. code-block:: Python import time import matplotlib.pyplot as plt import numpy as np from scipy.special import jn_zeros import magnelio as mio from magnelio import geo, ports from magnelio.constants import C0 .. GENERATED FROM PYTHON SOURCE LINES 35-45 Block 1 — before: one number drives both mesh scales ---------------------------------------------------- The mesh has two knobs: ``min_nodes_per_wavelength`` sets the bulk cell, ``min_cells_per_feature`` the cell at material interfaces. A ladder must scale *both*, in a fixed ratio — a ladder over the wavelength knob alone produces the same grid rung after rung as long as the interface cells dominate it, then jumps. Any other mesh setting your model needs (``min_cell_size``, ``max_edge_refinement``, …) goes in here too and stays fixed. .. GENERATED FROM PYTHON SOURCE LINES 45-57 .. code-block:: Python def rung(mnpw): """MeshControl for one rung: four bulk cells per interface cell.""" return mio.MeshControl( min_nodes_per_wavelength=mnpw, min_cells_per_feature=max(2, mnpw // 4), ) LADDER = (8, 12, 16, 24, 32, 48) .. GENERATED FROM PYTHON SOURCE LINES 58-69 Your simulation — wrapped in a function --------------------------------------- Whatever the script normally does between building the geometry and reading the result goes into one function that takes the ``MeshControl`` and returns the mesh and the result object. The example is a 50 Ω microstrip (the cross-section of tutorial 09, half-modelled on its symmetry plane) with a rectangular patch widening the trace in the middle: a shunt capacitance that reflects a good deal more than the −30 dB termination floor of a quasi-TEM port, so the S-parameters are the device's, not the port's. .. GENERATED FROM PYTHON SOURCE LINES 69-97 .. code-block:: Python h_sub, w_strip, t_strip = 0.8e-3, 1.2e-3, 0.2e-3 W_box, H_box, L = 8.0e-3, 5.0e-3, 16.0e-3 w_patch, l_patch = 4.0e-3, 3.0e-3 F_MAX = 15.0e9 fr4 = mio.Material.from_isotropic(name="FR4", epsilon=4.3) def simulate(mesh_control): substrate = geo.Brick(origin=(-W_box / 2, 0, 0), size=(W_box, h_sub, L), material=fr4) air = geo.Brick(origin=(-W_box / 2, h_sub, 0), size=(W_box, H_box - h_sub, L), material="air") strip = geo.Brick(origin=(-w_strip / 2, h_sub, 0), size=(w_strip, t_strip, L), material="pec") patch = geo.Brick( origin=(-w_patch / 2, h_sub, (L - l_patch) / 2), size=(w_patch, t_strip, l_patch), material="pec", ) trace = strip + patch model = mio.GeometryModel(boundary_conditions={"xmin": "SymmetryPMC"}) model.add(substrate) model.add(air - trace) model.add(trace) model.add_port(ports.PortWaveguide(name="port1", plane="zmin", n_modes=1)) model.add_port(ports.PortWaveguide(name="port2", plane="zmax", n_modes=1)) mesh = mio.Mesh.from_geometry(model, mesh_control, f_max=F_MAX) return mesh, mio.AnalysisScatteringTD(mesh=mesh, verbose=False).run(excited=["port1"]) .. GENERATED FROM PYTHON SOURCE LINES 98-120 Block 2 — after: the number and the stop rule --------------------------------------------- Every rung has its own frequency axis, so the S-parameters are first put on one common axis. The number compared is the largest magnitude of the *complex* difference between two consecutive rungs, over the band and over all S-parameters you care about: .. math:: \Delta S = \max_{f,\,ij} \; \bigl| S_{ij}^{(k)}(f) - S_{ij}^{(k-1)}(f) \bigr| and the stop rule is: :math:`\Delta S` below a threshold on two consecutive rungs. 0.02 is the customary threshold. It is the same scale a bench measurement lives on — 0.02 is a return-loss uncertainty of −34 dB, or 0.17 dB on a transmission near 0 dB — and, unlike a difference in decibels, it neither blows up at a transmission zero nor vanishes where a curve happens to be flat. The *first* rung of the passing pair is your production mesh: it is within tolerance of the next finer one, and the rung after that confirms the agreement was not a coincidence. ``converged_rung`` below is that rule. .. GENERATED FROM PYTHON SOURCE LINES 120-166 .. code-block:: Python F_COMMON = np.linspace(1.0e9, F_MAX, 141) PAIRS = [("port1", "port1"), ("port2", "port1")] TOL = 0.02 def s_table(result): """Complex S-parameters of `PAIRS` on the common axis, shape (pairs, f).""" return np.array([np.interp(F_COMMON, result.f_axis, result.S(*p)) for p in PAIRS]) def converged_rung(changes, tol): """Index of the first rung whose change *and* the next rung's are below tol. ``changes[k]`` is the change from rung k-1 to rung k (``nan`` for the first rung). Returns ``None`` when the ladder never settles. """ ok = [c < tol for c in changes] for k in range(len(ok) - 1): if ok[k] and ok[k + 1]: return k return None history = [] # (mnpw, S-table, ΔS to the previous rung) print("mnpw cells time max |ΔS| max Δ|S|") for mnpw in LADDER: t = time.perf_counter() mesh, result = simulate(rung(mnpw)) S = s_table(result) cells = mesh.Nx * mesh.Ny * mesh.Nz line = f"{mnpw:4d} {cells:6d} {time.perf_counter() - t:5.1f} s" delta = np.nan if history: delta = np.abs(S - history[-1][1]).max() delta_mag = np.abs(np.abs(S) - np.abs(history[-1][1])).max() line += f" {delta:8.4f} {delta_mag:8.4f} {'ok' if delta < TOL else ''}" print(line) history.append((mnpw, S, delta)) k = converged_rung([h[2] for h in history], TOL) if k is None: print(f"not converged within the ladder (ΔS < {TOL} on two consecutive rungs)") else: print(f"converged: mnpw {history[k][0]} (ΔS < {TOL} on this rung and the next)") .. rst-class:: sphx-glr-script-out .. code-block:: none mnpw cells time max |ΔS| max Δ|S| mesh | feature planes mesh | grid lines mesh | materials mesh | conformal cells mesh | PEC masks mesh | 11 x 17 x 15 cells 8 2805 0.7 s mesh | feature planes mesh | grid lines mesh | materials mesh | conformal cells mesh | PEC masks mesh | 15 x 21 x 22 cells 12 6930 1.3 s 0.1289 0.0216 mesh | feature planes mesh | grid lines mesh | materials mesh | conformal cells mesh | PEC masks mesh | 16 x 24 x 27 cells 16 10368 1.9 s 0.0431 0.0114 mesh | feature planes mesh | grid lines mesh | materials mesh | conformal cells mesh | PEC masks mesh | 21 x 27 x 42 cells 24 23814 3.9 s 0.0458 0.0129 mesh | feature planes mesh | grid lines mesh | materials mesh | conformal cells mesh | PEC masks mesh | 25 x 31 x 54 cells 32 41850 7.4 s 0.0158 0.0036 ok mesh | feature planes mesh | grid lines mesh | materials mesh | conformal cells mesh | PEC masks mesh | 32 x 38 x 81 cells 48 98496 27.1 s 0.0164 0.0050 ok converged: mnpw 32 (ΔS < 0.02 on this rung and the next) .. GENERATED FROM PYTHON SOURCE LINES 167-174 The curves behind the number ---------------------------- Top: the return loss of every rung on one axis — the family closes up as the ladder climbs. Bottom: :math:`|\Delta S|` against frequency for each pair of rungs, with the threshold drawn in. Where the deviation sits tells you what still moves. .. GENERATED FROM PYTHON SOURCE LINES 174-192 .. code-block:: Python fig, (ax_s, ax_d) = plt.subplots(2, 1, figsize=(6.4, 6.6), sharex=True) for mnpw, S, _ in history: ax_s.plot(F_COMMON / 1e9, 20 * np.log10(np.abs(S[0])), label=f"mnpw {mnpw}") ax_s.set_ylabel("|S11| [dB]") ax_s.set_title("capacitive patch on a microstrip, across the ladder") ax_s.grid(True, alpha=0.3) ax_s.legend(fontsize=8) for (m0, S0, _), (m1, S1, _) in zip(history, history[1:]): ax_d.plot(F_COMMON / 1e9, np.abs(S1 - S0).max(axis=0), label=f"{m0} → {m1}") ax_d.axhline(TOL, color="k", lw=0.8, ls="--", label=f"threshold {TOL}") ax_d.set_yscale("log") ax_d.set_xlabel("frequency [GHz]") ax_d.set_ylabel("max |ΔS| over S11, S21") ax_d.grid(True, alpha=0.3, which="both") ax_d.legend(fontsize=8, ncol=2) fig.tight_layout() .. image-sg:: /howto/images/sphx_glr_plot_mesh_convergence_001.png :alt: capacitive patch on a microstrip, across the ladder :srcset: /howto/images/sphx_glr_plot_mesh_convergence_001.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 193-221 What the number sees -------------------- Two things to know before trusting it. The deviation grows towards the top of the band and lives in the *phase* of S21: the magnitude-only column of the table is several times smaller than the complex one. That is the numerical dispersion of the feed lines — 16 mm of line is one and a half wavelengths at 15 GHz, and the grid propagates a wave a little too slowly on a coarse mesh. A ladder on the bare line without the patch gives nearly the same :math:`\Delta S` as this one. The complex criterion is right to count phase (a filter or a coupler depends on it), but it means that long leads to the ports set the rung, not the device. Keep the leads as short as the ports allow. The ladder is not monotone: the deviation went *up* from the 12 → 16 pair to the 16 → 24 pair. On a conformal grid the error depends on where the grid lines happen to cut the trace edges and the substrate, and that changes from rung to rung. This is why the stop rule asks for two consecutive passes and not for one. And a caveat that belongs to strongly resonant structures: a small shift of a resonance moves the S-parameters at that frequency by a lot, so :math:`\Delta S` stays large long after the *resonance frequency* has converged to a fraction of a percent. A resonator or a narrow-band filter is judged on its frequencies — the second half of this page. .. GENERATED FROM PYTHON SOURCE LINES 223-235 The same ladder around a resonance ---------------------------------- Blocks 1 and 2 again, with the number now a frequency and the stop rule its relative change. The model is a cylindrical metal cavity (pillbox), because its lowest mode has a closed form to hold the ladder against: :math:`f_{010} = 2.4048\,c_0 / (2\pi R)`, independent of the height. The curved wall is where the grid has to work — it is represented by the conformal sub-cell treatment of partially filled cells, not by a staircase. The ladder starts at 12 because a driven run tolerates a coarser grid than an eigenmode of a curved cavity. .. GENERATED FROM PYTHON SOURCE LINES 235-277 .. code-block:: Python R_CAV, H_CAV = 10e-3, 6e-3 F_EXACT = float(jn_zeros(0, 1)[0]) * C0 / (2 * np.pi * R_CAV) TOL_F = 1e-2 # relative change in f0 that your specification allows def simulate_cavity(mesh_control): model = mio.GeometryModel(background="pec") model.add(geo.Cylinder(origin=(0, 0, 0), radius=R_CAV, height=H_CAV, axis="z", material="air")) mesh = mio.Mesh.from_geometry(model, mesh_control, f_max=14e9) result = mio.AnalysisEigenmode( mesh=mesh, n_modes=2, sigma=(2 * np.pi * 11e9) ** 2, verbose=False ).run() return float(result.frequencies[0]), mesh.Nx * mesh.Ny * mesh.Nz ladder_f = (12, 16, 24, 32, 48) f0, changes = [], [] print(f"exact TM010: {F_EXACT / 1e9:.4f} GHz") print("mnpw f0 [GHz] cells time change error") for mnpw in ladder_f: t = time.perf_counter() f, cells = simulate_cavity(rung(mnpw)) change = np.nan if not f0 else abs(f - f0[-1]) / f0[-1] print( f"{mnpw:4d} {f / 1e9:8.4f} {cells:6d} {time.perf_counter() - t:4.1f} s " f"{'' if not f0 else f'{change * 100:.2f} %':8s} " f"{(f - F_EXACT) / F_EXACT * 100:+.2f} % {'ok' if change < TOL_F else ''}" ) f0.append(f) changes.append(change) f0 = np.array(f0) k = converged_rung(changes, TOL_F) if k is None: print(f"not converged within the ladder ({TOL_F * 100:.1f} % on two consecutive rungs)") else: print( f"converged: mnpw {ladder_f[k]} at {TOL_F * 100:.1f} % — " f"its actual error is {abs(f0[k] - F_EXACT) / F_EXACT * 100:.2f} %" ) .. rst-class:: sphx-glr-script-out .. code-block:: none exact TM010: 11.4743 GHz mnpw f0 [GHz] cells time change error mesh | feature planes mesh | grid lines mesh | materials mesh | conformal cells mesh | PEC masks mesh | 12 x 12 x 4 cells 12 11.4057 576 0.1 s -0.60 % mesh | feature planes mesh | grid lines mesh | materials mesh | conformal cells mesh | PEC masks mesh | 15 x 15 x 5 cells 16 11.2897 1125 0.1 s 1.02 % -1.61 % mesh | feature planes mesh | grid lines mesh | materials mesh | conformal cells mesh | PEC masks mesh | 23 x 23 x 7 cells 24 11.3901 3703 0.4 s 0.89 % -0.73 % ok mesh | feature planes mesh | grid lines mesh | materials mesh | conformal cells mesh | PEC masks mesh | 30 x 30 x 9 cells 32 11.4487 8100 2.2 s 0.51 % -0.22 % ok mesh | feature planes mesh | grid lines mesh | materials mesh | conformal cells mesh | PEC masks mesh | 45 x 45 x 14 cells 48 11.4889 28350 40.4 s 0.35 % +0.13 % ok converged: mnpw 24 at 1.0 % — its actual error is 0.73 % .. GENERATED FROM PYTHON SOURCE LINES 278-298 Reading the ladder — and extrapolating it ----------------------------------------- The *change between rungs* against your tolerance is the criterion, as before — and here the exact value checks the rule: at a 1 % tolerance it settles on rung 24, whose actual error is 0.7 %. A resonance ladder offers one thing more: a second-order scheme converges like :math:`h^2`, so plotting the frequency against :math:`1/N^2` (:math:`N` the nodes per wavelength) gives a straight line once the mesh is in the asymptotic regime, and the intercept at :math:`1/N^2 \to 0` is an estimate of the converged value. Here the exact value is known, so the plot also shows how good that estimate is — and it is a caution as much as a demonstration: the coarse rungs wobble with the way the grid lines cut the curved wall, the fine rungs approach from below and slightly overshoot, and the extrapolation lands a few tenths of a percent off. A conformal boundary carries an error component that does not scale as :math:`h^2`. Take the extrapolation as an estimate and the change between rungs as the criterion. .. GENERATED FROM PYTHON SOURCE LINES 298-318 .. code-block:: Python x = 1.0 / np.asarray(ladder_f, dtype=float) ** 2 slope, intercept = np.polyfit(x[-3:], f0[-3:], 1) # fit the three finest rungs fig, ax = plt.subplots(figsize=(6.4, 4.2)) ax.plot(x, f0 / 1e9, "o-", label="ladder") xx = np.linspace(0.0, x.max(), 50) ax.plot(xx, (slope * xx + intercept) / 1e9, "--", label="fit through the three finest rungs") ax.plot(0.0, intercept / 1e9, "s", label=f"extrapolated {intercept / 1e9:.4f} GHz") ax.axhline(F_EXACT / 1e9, color="k", lw=0.8, label=f"exact {F_EXACT / 1e9:.4f} GHz") ax.set_xlabel("1 / (nodes per wavelength)²") ax.set_ylabel("f₀ [GHz]") ax.set_title("TM₀₁₀ of a pillbox versus resolution") ax.grid(True, alpha=0.3) ax.legend() print( f"extrapolated f0: {intercept / 1e9:.4f} GHz ({(intercept - F_EXACT) / F_EXACT * 1e6:+.0f} ppm)" ) .. image-sg:: /howto/images/sphx_glr_plot_mesh_convergence_002.png :alt: TM₀₁₀ of a pillbox versus resolution :srcset: /howto/images/sphx_glr_plot_mesh_convergence_002.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none extrapolated f0: 11.5223 GHz (+4184 ppm) .. GENERATED FROM PYTHON SOURCE LINES 319-330 Two more things a resonance ladder teaches. A high-contrast dielectric body (a ceramic puck with εᵣ = 45) run through the same rungs moves by ±3 % *in either direction* between them — no line to extrapolate, and the spread across the rungs is the honest error of its absolute frequency. And a *ratio* of two results from the same mesh — a coupling coefficient from two eigenfrequencies, a Q, a relative bandwidth — converges much faster than either of them, because both carry nearly the same discretisation error. Tutorial 13 designs a filter on that basis. Run the ladder on the quantity you actually design with, and read it against that quantity's own tolerance. .. GENERATED FROM PYTHON SOURCE LINES 332-352 What the ladder cannot see -------------------------- A ladder finds errors that *shrink* with the cell size. It is blind to a feature the grid does not contain at all: a chamfer or fillet smaller than half a cell has no effect on the result, then appears in one step when a rung finally resolves it — the ladder shows a jump, not a trend. The mesher places grid planes on such edges so the feature occupies a cell layer of its own, and warns when an edge would need a cell finer than ``h_max / max_edge_refinement`` and is dropped. Read those warnings before trusting a ladder; the remedy they name (``MeshControl(max_edge_refinement=...)``) belongs into ``rung()`` as a fixed setting, not into the resolution. The flip side: a feature layer the mesher *does* resolve keeps its size — one cell across the chamfer, whatever the rung — so the grid is no longer refined uniformly and the :math:`h^2` trend is disturbed. Run the ladder on the sharp-edged model to judge the resolution, then add the features back. The mechanism is explained under :doc:`../methods/meshing-conformal`. .. GENERATED FROM PYTHON SOURCE LINES 352-354 .. code-block:: Python plt.show() .. rst-class:: sphx-glr-timing **Total running time of the script:** (1 minutes 25.988 seconds) .. _sphx_glr_download_howto_plot_mesh_convergence.py: .. only:: html .. container:: sphx-glr-footer sphx-glr-footer-example .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: plot_mesh_convergence.ipynb ` .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: plot_mesh_convergence.py ` .. container:: sphx-glr-download sphx-glr-download-zip :download:`Download zipped: plot_mesh_convergence.zip ` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_