# FEM Dataset Lookup in SAX Models Circuit models of superconducting layouts need electrostatic quantities such as capacitance matrices, which for realistic geometries come from finite-element (FEM) solvers [[Jin14](#id34)]. One solve takes minutes to hours, far too slow to call inside a circuit simulation or an optimisation loop. So a sweep over the design parameters is solved once, stored as a dataset, and the circuit model interpolates between the stored points. [`qpdk.models.datasets`](https://gdsfactory.github.io/quantum-rf-pdk/simulation.html.md#module-qpdk.models.datasets) stores the results as a long-format [Polars](https://pola.rs) table with one row per matrix entry. Each Parquet file (or the Delta table schema) also carries the dataset metadata: units, terminal and ground conventions, whether the data is synthetic, and free-form provenance such as the layer stack and solver. A file is therefore self-describing, and the parameter grid is read off the data rather than declared separately. Curated datasets ship in the package as Parquet parts in Git LFS; shared, growing datasets can live in a Delta Lake table on a cloud bucket. This notebook follows generation, storage, lazy queries and JAX interpolation, then uses the capacitance lookup in a SAX model. Run `uv run --script notebooks/src/fem_dataset_lookup.py` with Git LFS installed to use the pinned QPDK revision and fetch its bundled results automatically. The bundled data comes from Palace electrostatic solves of the actual QPDK metal polygons, including their short leads, as perfect conductor sheets on silicon. A separate coplanar ground frame surrounds the device. The geometry and mesh settings travel with the data. It covers lengths 40–120 µm, widths 5–20 µm and gaps 4–10 µm; these results apply to that geometry and stack. ```ipython3 import sys if "google.colab" in sys.modules: import subprocess print("Running in Google Colab. Installing QPDK...") subprocess.check_call(["apt-get", "install", "-y", "-qq", "git-lfs"]) subprocess.check_call(["git", "lfs", "install", "--skip-repo"]) subprocess.check_call([ sys.executable, "-m", "pip", "install", "-q", "qpdk[models] @ git+https://github.com/jackgdsf/quantum-rf-pdk.git@b9acce2774c9cdbcc20ad280916ed76265a7a74a", "matplotlib-inline", ]) ``` ```ipython3 import subprocess import time from pathlib import Path import jax import jax.numpy as jnp import matplotlib.pyplot as plt import numpy as np import polars as pl import sax import xarray as xr from matplotlib.colors import LogNorm from matplotlib_inline.backend_inline import set_matplotlib_formats from qpdk import logger from qpdk.colormaps import QPDK from qpdk.config import PATH from qpdk.models.capacitor import plate_capacitor from qpdk.models.constants import DEFAULT_FREQUENCY from qpdk.models.couplers import ( cpw_coupling_model, cpw_cpw_coupling_capacitance, cpw_cpw_coupling_capacitance_per_length_analytical, ) from qpdk.models.cpw import cpw_z0_from_cross_section from qpdk.models.datasets import ( Dataset, GridInterpolator, capacitance_model, check_maxwell, maxwell_to_mutual, ) from qpdk.models.datasets.generate import write from qpdk.models.generic import capacitor from qpdk.tech import coplanar_waveguide set_matplotlib_formats("png", "svg") for style_source in (PATH.repo / "docs" / "qpdk.mplstyle", "qpdk"): try: plt.style.use(style_source) except OSError: continue break ``` ## Inspect the dataset The metadata is read from the files and validated when the dataset opens. The results table is a plain Polars frame with one row per matrix entry, so it can be filtered and aggregated without loading anything into JAX. ```ipython3 dataset = Dataset("plate_capacitor_palace") metadata = dataset.metadata logger.info(f"{dataset!r}: synthetic={metadata.synthetic}") logger.info(f"solver: {metadata.provenance['solver']}") logger.info(f"extraction settings: {metadata.provenance['settings']}") logger.info(f"terminals {metadata.terminals} w.r.t. {metadata.reference_ground!r}") for axis in metadata.axes: logger.info( f"{axis.name} [{axis.unit}]: validated domain {axis.validated or 'grid span'}" ) ``` ```myst-ansi 2026-10-22 13:29:39.435 | INFO  | __main__::3 - Dataset('plate_capacitor_palace', store=ParquetParts('/home/runner/work/quantum-rf-pdk/quantum-rf-pdk/qpdk/models/datasets/data/plate_capacitor_palace')): synthetic=False ``` ```myst-ansi 2026-10-22 13:29:39.436 | INFO  | __main__::4 - solver: Palace ``` ```myst-ansi 2026-10-22 13:29:39.437 | INFO  | __main__::5 - extraction settings: {'near_mesh': 0.55, 'far_mesh': 25.0, 'domain_pad': 100.0, 'ground_clearance': 10.0, 'ground_pad': 45.0, 'permittivity': 11.45, 'order': 2, 'tolerance': 1e-09, 'save_fields': False} ``` ```myst-ansi 2026-10-22 13:29:39.437 | INFO  | __main__::6 - terminals ('o1', 'o2') w.r.t. 'coplanar_ground_frame' ``` ```myst-ansi 2026-10-22 13:29:39.438 | INFO  | __main__::8 - length [um]: validated domain grid span ``` ```myst-ansi 2026-10-22 13:29:39.438 | INFO  | __main__::8 - width [um]: validated domain grid span ``` ```myst-ansi 2026-10-22 13:29:39.439 | INFO  | __main__::8 - gap [um]: validated domain grid span ``` ```ipython3 dataset.table.head(8) ``` ```myst-ansi shape: (8, 12) ┌──────────────────┬────────┬────────┬───────┬───┬──────┬─────────────┬────────────┬──────┐ │ run_id ┆ status ┆ length ┆ width ┆ … ┆ col ┆ value ┆ value_imag ┆ unit │ │ --- ┆ --- ┆ --- ┆ --- ┆ ┆ --- ┆ --- ┆ --- ┆ --- │ │ str ┆ str ┆ f64 ┆ f64 ┆ ┆ str ┆ f64 ┆ f64 ┆ str │ ╞══════════════════╪════════╪════════╪═══════╪═══╪══════╪═════════════╪════════════╪══════╡ │ b42beed1e26f12ae ┆ ok ┆ 40.0 ┆ 5.0 ┆ … ┆ null ┆ 0.281896 ┆ null ┆ 1 │ │ b42beed1e26f12ae ┆ ok ┆ 40.0 ┆ 5.0 ┆ … ┆ o1 ┆ 6.5721e-15 ┆ null ┆ F │ │ b42beed1e26f12ae ┆ ok ┆ 40.0 ┆ 5.0 ┆ … ┆ o2 ┆ -2.2157e-15 ┆ null ┆ F │ │ b42beed1e26f12ae ┆ ok ┆ 40.0 ┆ 5.0 ┆ … ┆ o1 ┆ -2.2157e-15 ┆ null ┆ F │ │ b42beed1e26f12ae ┆ ok ┆ 40.0 ┆ 5.0 ┆ … ┆ o2 ┆ 6.5759e-15 ┆ null ┆ F │ │ b42beed1e26f12ae ┆ ok ┆ 40.0 ┆ 5.0 ┆ … ┆ null ┆ 12.0 ┆ null ┆ 1 │ │ b42beed1e26f12ae ┆ ok ┆ 40.0 ┆ 5.0 ┆ … ┆ null ┆ 3.2135e-10 ┆ null ┆ 1 │ │ a2293eecb392f6e3 ┆ ok ┆ 40.0 ┆ 5.0 ┆ … ┆ null ┆ 0.271239 ┆ null ┆ 1 │ └──────────────────┴────────┴────────┴───────┴───┴──────┴─────────────┴────────────┴──────┘ ``` ```ipython3 ( dataset .scan() .filter(pl.col("row") != pl.col("col")) .group_by("gap") .agg(c_mutual_fF=(-pl.col("value")).mean() * 1e15) .sort("gap") .collect() ) ``` ```myst-ansi shape: (3, 2) ┌──────┬─────────────┐ │ gap ┆ c_mutual_fF │ │ --- ┆ --- │ │ f64 ┆ f64 │ ╞══════╪═════════════╡ │ 4.0 ┆ 5.735053 │ │ 7.0 ┆ 4.348927 │ │ 10.0 ┆ 3.538881 │ └──────┴─────────────┘ ``` ## Define a generation script `qpdk/models/datasets/data/plate_capacitor.py` is a complete experiment: its `GRID`, `SETTINGS`, metadata and `solve` function describe what to extract. Its inline dependencies let `uv run --script` create the Python 3.12 environment required by gsim. Copy the script beside the original in `qpdk/models/datasets/data/` and edit it for another dataset. Meshing runs on the host, including when Palace uses `--sif`. Gmsh needs its native libraries there; on Linux this includes `libGLU.so.1`. The container needs a Palace binary with MPI support; select its executable with `--container-binary`. Preview the full sweep or run it from a checkout: ```bash uv run --script qpdk/models/datasets/data/plate_capacitor.py --dry-run uv run --script qpdk/models/datasets/data/plate_capacitor.py --processes 4 # A local Palace container: uv run --script qpdk/models/datasets/data/plate_capacitor.py --sif /path/to/palace.sif ``` `--processes` sets CPU threads for meshing and MPI ranks for Palace. The script builds a QPDK component. gsim meshes its electrodes with meshwell, writes the configuration and runs Palace. Results carry their metadata; lookup needs neither the generator nor the simulator. Generator scripts stay outside wheels. ```ipython3 small_grid = {"length": [40.0, 80.0], "width": [10.0], "gap": [4.0, 7.0, 10.0]} output = Path("build/dataset-example/coarse") generator_script = Path("qpdk/models/datasets/data/plate_capacitor.py") ``` ## Run a small sweep and inspect its output Copy the generator beside the original, replace `GRID` with `small_grid`, then run it with `--output build/dataset-example/coarse`. For example, its grid definition becomes: ```python GRID = {"length": [40.0, 80.0], "width": [10.0], "gap": [4.0, 7.0, 10.0]} ``` Set `RUN_PALACE=True` to execute the generator from a checkout. It runs the grid defined in `generator_script`; point this path to your edited copy. The default selects six existing Palace results and performs no new solves. Both paths produce a self-describing table with four matrix entries and scalar diagnostics per solve. ```ipython3 RUN_PALACE = False if RUN_PALACE: subprocess.run( # ruff: ignore[subprocess-without-shell-equals-true] [ "uv", "run", "--script", str(generator_script), "--output", str(output), ], check=True, ) generated = Dataset(output) else: selection = ( dataset .scan() .filter(*(pl.col(axis).is_in(values) for axis, values in small_grid.items())) .collect() ) generated = write(output, dataset.metadata, selection) logger.info("Selected six existing Palace solves; no simulator was run") coarse = generated.grid("maxwell_capacitance", cross_section="cpw") logger.info(f"Output: {generated.scan().select(pl.len()).collect().item()} rows") logger.info( f"Dense grid: {coarse.values.shape}, terminals={generated.metadata.terminals}" ) coarse_lookup = GridInterpolator(coarse) logger.info( f"Small-sweep lookup at 60/10/7 µm [fF]:\n" f"{coarse_lookup(length=60.0, width=10.0, gap=7.0) * 1e15}" ) ``` ```myst-ansi 2026-10-22 13:29:39.504 | INFO  | __main__::24 - Selected six existing Palace solves; no simulator was run ``` ```myst-ansi 2026-10-22 13:29:39.517 | INFO  | __main__::26 - Output: 42 rows ``` ```myst-ansi 2026-10-22 13:29:39.517 | INFO  | __main__::27 - Dense grid: (2, 1, 3, 2, 2), terminals=('o1', 'o2') ``` ```myst-ansi 2026-10-22 13:29:40.576 | INFO  | __main__::31 - Small-sweep lookup at 60/10/7 µm [fF]: [[10.89142216 -3.13709128] [-3.13709128 10.90180517]] ``` ```ipython3 ( generated .scan() .filter(pl.col("quantity") == "maxwell_capacitance") .group_by("length", "width", "gap", "status") .len(name="matrix_entries") .sort("length", "gap") .collect() ) ``` ```myst-ansi shape: (6, 5) ┌────────┬───────┬──────┬────────┬────────────────┐ │ length ┆ width ┆ gap ┆ status ┆ matrix_entries │ │ --- ┆ --- ┆ --- ┆ --- ┆ --- │ │ f64 ┆ f64 ┆ f64 ┆ str ┆ u32 │ ╞════════╪═══════╪══════╪════════╪════════════════╡ │ 40.0 ┆ 10.0 ┆ 4.0 ┆ ok ┆ 4 │ │ 40.0 ┆ 10.0 ┆ 7.0 ┆ ok ┆ 4 │ │ 40.0 ┆ 10.0 ┆ 10.0 ┆ ok ┆ 4 │ │ 80.0 ┆ 10.0 ┆ 4.0 ┆ ok ┆ 4 │ │ 80.0 ┆ 10.0 ┆ 7.0 ┆ ok ┆ 4 │ │ 80.0 ┆ 10.0 ┆ 10.0 ┆ ok ┆ 4 │ └────────┴───────┴──────┴────────┴────────────────┘ ``` ## Resume and extend a sweep Rerun the generator after an interruption. Each geometry retains its inputs, mesh, config, solver log and output. Completed matching runs are reused; failed ones retry. Geometry, mesh, generator or runtime changes get a new run directory. The published dataset is replaced only after the sweep succeeds. Add `120.0` to the small script’s length grid and rerun with the same `--workdir`: six geometries are reused and three are added. Use `--output` and `--workdir` to select locations; paths are relative to your current directory. ## Parallel sweeps Both experiments accept `--shard INDEX --shards COUNT`. The longest axis is split into disjoint rectangular grids, with unchanged solver inputs and run identities. Workers write separate outputs and run directories; merging checks provenance, coverage and successful results before publication. Use a small array of longer shards, with no more tasks than points on the longest axis. For the CPW grid: ```bash generator=qpdk/models/datasets/data/cpw_coupling.py sbatch --array=0-1 --cpus-per-task=4 qpdk/models/datasets/data/slurm_array.sh "$generator" \ build/cpw-shards build/cpw-runs --sif /path/to/palace.sif # After every task succeeds, pass each shard as a repeated option: shards=() for i in {0..1}; do shards+=(--merge-shards "build/cpw-shards/shard-$i"); done uv run --script "$generator" "${shards[@]}" --output build/datasets/cpw_coupling_palace ``` Each task runs one solve at a time with four MPI ranks, one thread per rank. Specify the array size with `--array`; `--cpus-per-task` sets MPI ranks per solve. Without Slurm, run the entire script directly or use a local loop. This sequential loop needs four cores; independent shard commands can also run concurrently if each has its own four-core allocation: ```bash generator=qpdk/models/datasets/data/cpw_coupling.py export OMP_NUM_THREADS=1 OPENBLAS_NUM_THREADS=1 MKL_NUM_THREADS=1 shards=() for i in {0..1}; do uv run --script "$generator" --shard "$i" --shards 2 --processes 4 \ --output "build/cpw-shards/shard-$i" --workdir "build/cpw-runs/shard-$i" shards+=(--merge-shards "build/cpw-shards/shard-$i") done uv run --script "$generator" "${shards[@]}" --output build/datasets/cpw_coupling_palace ``` Keep the script, settings, runtime and array size fixed while resuming. The reusable `partition_grid` and `merge` functions also work with other experiment scripts and schedulers. ## Check mesh and domain sensitivity Each result stores the worst terminal’s linear-solver residual and iteration count, plus Palace’s energy-normalized recovered-flux FEM error indicator. The generator rejects missing diagnostics or a linear solve that exceeds `SETTINGS.tolerance`. Linear convergence alone does not establish mesh accuracy. The FEM indicator ranks candidates for refinement; it is not a percentage error in capacitance. A denser parameter grid cannot repair a coarse FEM mesh. In a copy of the generator, select one geometry and vary the settings: ```python from dataclasses import replace GRID = {"length": [80.0], "width": [10.0], "gap": [7.0]} SETTINGS = replace(SETTINGS, near_mesh=0.38, far_mesh=17.0, save_fields=True) ``` Reduce both near and far mesh sizes with separate output directories, then increase `domain_pad` to 150 µm. Compare the Maxwell matrices and inspect the saved fields in ParaView, especially near metal edges and gaps. Repeat representative checks near the ends of any expanded geometry range. For each quantity of interest, compare successive solves using $|C_{\text{fine}}-C_{\text{coarse}}|/|C_{\text{fine}}|$, checking mesh and domain changes separately. The bundled plate dataset records representative finer-mesh and larger-domain comparisons, a native-runtime cross-check and fresh off-grid references in its provenance. These finite checks do not establish a uniform accuracy bound. ```ipython3 plate_checks = metadata.provenance["validation"] logger.info( f"Native runtime grid comparison: {plate_checks['native_grid_max_relative']:.3%}; " f"{plate_checks['held_out_points']} off-grid references: " f"maximum lookup difference {plate_checks['held_out_max_relative']:.2%}" ) ( pl .DataFrame(plate_checks["mesh_domain"]) .unnest("point") .select( "length", "width", "gap", (100 * pl.col("refined_relative")).alias("Finer mesh (%)"), (100 * pl.col("domain_relative")).alias("Larger domain (%)"), ) ) ``` ```myst-ansi 2026-10-22 13:29:40.592 | INFO  | __main__::2 - Native runtime grid comparison: 0.000%; 4 off-grid references: maximum lookup difference 2.79% ``` ```myst-ansi shape: (3, 5) ┌────────┬───────┬──────┬────────────────┬───────────────────┐ │ length ┆ width ┆ gap ┆ Finer mesh (%) ┆ Larger domain (%) │ │ --- ┆ --- ┆ --- ┆ --- ┆ --- │ │ f64 ┆ f64 ┆ f64 ┆ f64 ┆ f64 │ ╞════════╪═══════╪══════╪════════════════╪═══════════════════╡ │ 40.0 ┆ 5.0 ┆ 4.0 ┆ 0.73666 ┆ 0.1187 │ │ 80.0 ┆ 10.0 ┆ 7.0 ┆ 0.782793 ┆ 0.168654 │ │ 120.0 ┆ 20.0 ┆ 10.0 ┆ 0.776481 ┆ 0.282732 │ └────────┴───────┴──────┴────────────────┴───────────────────┘ ``` ## Query a shared Delta table Cloud datasets use the same interface. Select the quantity and variant before collecting; `grid()` does this internally. Pin a version for reproducible models: ```python shared = Dataset("gs://my-bucket/capacitors", delta=True, version=3) rows = shared.scan().filter( pl.col("quantity") == "maxwell_capacitance", pl.col("cross_section") == "cpw", ).select("length", "gap", "value").collect() lookup = GridInterpolator(shared.grid("maxwell_capacitance", cross_section="cpw")) ``` Scans stay lazy until `collect()`, allowing Polars to push filters and column selection into storage. Only the selected interpolation grid becomes a dense in-memory array. `dataset.table` explicitly loads and validates all rows. ## Look up a Maxwell capacitance matrix For conductors $1, \dots, n$ at potentials $V_j$ above a reference ground, the charges are linear in the potentials, $Q_i = \sum_j C_{ij} V_j$. The *Maxwell* capacitance matrix $C$ is what an electrostatic FEM solver returns: it is symmetric, its diagonal is positive, and its off-diagonal entries are non-positive. The *mutual* capacitance between conductors $i \neq j$ is $C^{\text{m}}_{ij} = -C_{ij}$, and the capacitance of conductor $i$ to ground is the row sum $C^{\text{m}}_{ii} = \sum_j C_{ij}$. The mutual form is what appears as branch capacitances in a lumped circuit, and the Maxwell form is the capacitance matrix of circuit quantization [[VD17](#id71)]. `maxwell_to_mutual()` converts between them, and `check_maxwell()` raises if any property above is violated. Between the solved points the lookup is multilinear: inside the grid cell containing a query, the result is a weighted mean of the $2^n$ cell corners, with weights that are products of the 1D linear weights along each axis [[WZ88](#id72)]. It reproduces the data exactly at grid points, is continuous everywhere, and needs no tuning. Outside the validated domain there is no data to weigh, so the result is NaN. ![Rectilinear grid of solved points; a query inside a cell is a weighted mean of the cell corners, a query outside the grid returns NaN](notebooks/figures/fem-dataset-lookup.svg) `grid()` arranges the table into a dense `(length, width, gap, 2, 2)` array and raises if a grid point is missing or failed. Discrete variants such as the cross-section are chosen explicitly here and never interpolated. `GridInterpolator` wraps [`jax.scipy.interpolate.RegularGridInterpolator`](https://docs.jax.dev/en/latest/_autosummary/jax.scipy.interpolate.RegularGridInterpolator.html#jax.scipy.interpolate.RegularGridInterpolator) and holds plain JAX arrays, so it works inside [`jax.jit()`](https://docs.jax.dev/en/latest/_autosummary/jax.jit.html#jax.jit), [`jax.vmap()`](https://docs.jax.dev/en/latest/_autosummary/jax.vmap.html#jax.vmap), and [`jax.grad()`](https://docs.jax.dev/en/latest/_autosummary/jax.grad.html#jax.grad). ```ipython3 grid = dataset.grid("maxwell_capacitance", cross_section="cpw") c_maxwell = GridInterpolator(grid) logger.info(repr(c_maxwell)) C = c_maxwell(length=100.0, width=10.0, gap=5.0) logger.info(f"Maxwell matrix [fF]:\n{C * 1e15}") logger.info(f"Mutual matrix [fF]:\n{maxwell_to_mutual(C) * 1e15}") check_maxwell(C) # raises NonPhysicalMatrixError if C is not physical ``` ```myst-ansi 2026-10-22 13:29:40.623 | INFO  | __main__::3 - GridInterpolator('plate_capacitor_palace', 'maxwell_capacitance', domain={'length': (40.0, 120.0), 'width': (5.0, 20.0), 'gap': (4.0, 10.0)}, out_of_range='nan') ``` ```myst-ansi 2026-10-22 13:29:40.751 | INFO  | __main__::6 - Maxwell matrix [fF]: [[17.73218645 -6.60501633] [-6.60501633 17.71842436]] ``` ```myst-ansi 2026-10-22 13:29:40.990 | INFO  | __main__::7 - Mutual matrix [fF]: [[11.12717012 6.60501633] [ 6.60501633 11.11340803]] ``` Outside the validated domain the lookup returns NaN. [`sax.interpolate_xarray()`](https://gdsfactory.github.io/sax/api/#sax.interpolate_xarray) clamps to the nearest edge instead, which hides an out-of-range design. ```ipython3 lengths = jnp.linspace(0.0, 250.0, 251) ours = -c_maxwell(length=lengths, width=10.0, gap=5.0)[:, 0, 1] xarr = xr.DataArray( grid.values.reshape(*grid.values.shape[:3], 4), coords={ **dict(zip(grid.axis_names, grid.coords, strict=True)), "targets": np.array(["c11", "c12", "c21", "c22"], dtype=object), }, ) theirs = -jax.vmap( lambda length: sax.interpolate_xarray( xarr, length=length, width=jnp.asarray(10.0), gap=jnp.asarray(5.0) )["c12"] )(lengths) fig, ax = plt.subplots() ax.plot(lengths, theirs * 1e15, "--", label="sax.interpolate_xarray (clamps)") ax.plot(lengths, ours * 1e15, label="GridInterpolator (NaN outside)") ax.axvspan(*c_maxwell.domain["length"], alpha=0.1, label="validated domain") ax.set_xlabel(r"Length $l$ ($\text{µm}$)") ax.set_ylabel(r"Mutual capacitance $C_{12}$ ($\text{fF}$)") ax.legend() plt.show() ``` ```myst-ansi
``` ## Replace an explicit capacitance in a SAX model `plate_capacitor()` computes its capacitance from a closed-form expression. The same circuit can take the mutual capacitance from the dataset instead. The analytical formula and FEM model include different geometry and ground assumptions, so their disagreement is visible here. The lookup uses only the mutual branch; a complete circuit can also include the extracted pad-to-ground capacitances. `capacitance_model` supplies that complete lumped N-port network for any number of terminals; `s_parameters_model` directly interpolates stored N-port scattering data. ```ipython3 def plate_capacitor_lookup( *, f=DEFAULT_FREQUENCY, length: float = 80.0, width: float = 5.0, gap: float = 7.0, cross_section="cpw", ) -> sax.SDict: """Plate capacitor with the mutual capacitance looked up from the dataset.""" f = jnp.asarray(f) c_mutual = maxwell_to_mutual(c_maxwell(length=length, width=width, gap=gap))[0, 1] return capacitor( f=f, capacitance=c_mutual, z0=cpw_z0_from_cross_section(cross_section, f) ) complete_plate = capacitance_model(dataset, cross_section="cpw") logger.info( f"Complete two-terminal network at 5 GHz: " f"{complete_plate(f=5e9, length=80.0, width=10.0, gap=7.0)}" ) f = jnp.linspace(1e9, 10e9, 201) s_lookup = jax.jit( lambda gap: plate_capacitor_lookup(f=f, length=80.0, width=10.0, gap=gap) )(7.0) s_analytical = plate_capacitor(f=f, length=80.0, width=10.0, gap=7.0) fig, ax = plt.subplots() ax.plot(f / 1e9, 20 * jnp.log10(jnp.abs(s_analytical["o1", "o2"])), label="analytical") ax.plot( f / 1e9, 20 * jnp.log10(jnp.abs(s_lookup["o1", "o2"])), "--", label="dataset lookup" ) ax.set_xlabel(r"Frequency $f$ ($\text{GHz}$)") ax.set_ylabel(r"$|S_{21}|$ ($\text{dB}$)") ax.legend() plt.show() ``` ```myst-ansi 2026-10-22 13:29:46.130 | INFO  | __main__::18 - Complete two-terminal network at 5 GHz: {('o1', 'o1'): Array(0.9989619-0.04352421j, dtype=complex128), ('o1', 'o2'): Array(0.00058555+0.01343248j, dtype=complex128), ('o2', 'o1'): Array(0.00058555+0.01343248j, dtype=complex128), ('o2', 'o2'): Array(0.99895995-0.04356874j, dtype=complex128)} ``` ```myst-ansi
``` The lookup is differentiable, so the gap can be tuned by gradient descent towards a target transmission. ```ipython3 d_s21_d_gap = jax.jit( jax.grad( lambda gap: jnp.abs( plate_capacitor_lookup(f=5e9, length=80.0, width=10.0, gap=gap)["o1", "o2"] ) ) ) logger.info(f"d|S21|/d gap at 5 GHz, gap = 7 µm: {float(d_s21_d_gap(7.0)):.3e} 1/µm") ``` ```myst-ansi 2026-10-22 13:29:48.445 | INFO  | __main__::8 - d|S21|/d gap at 5 GHz, gap = 7 µm: -1.448e-03 1/µm ``` ## Compile and query time Compile time and steady-state query time are measured separately, for a single point and for a batch of 10 000 points. ```ipython3 def benchmark(fn, *args, repeats: int = 50) -> tuple[float, float]: """Return (compile time, steady-state time per call) in ms.""" jitted = jax.jit(fn) t0 = time.perf_counter() jax.block_until_ready(jitted(*args)) compile_ms = (time.perf_counter() - t0) * 1e3 t0 = time.perf_counter() for _ in range(repeats): jax.block_until_ready(jitted(*args)) return compile_ms, (time.perf_counter() - t0) / repeats * 1e3 rng = np.random.default_rng(0) batch = { name: jnp.asarray(rng.uniform(lo, hi, 10_000)) for name, (lo, hi) in c_maxwell.domain.items() } single = {name: x[0] for name, x in batch.items()} rows = [] for label, point in [("1 point", single), ("10 000 points", batch)]: for method, fn in [ ("GridInterpolator", lambda p: c_maxwell(**p)[..., 0, 1]), ("sax.interpolate_xarray", lambda p: sax.interpolate_xarray(xarr, **p)["c12"]), ]: compile_ms, query_ms = benchmark(fn, point) rows.append({ "method": method, "query": label, "compile_ms": compile_ms, "query_ms": query_ms, }) pl.DataFrame(rows) ``` ```myst-ansi shape: (4, 4) ┌────────────────────────┬───────────────┬────────────┬──────────┐ │ method ┆ query ┆ compile_ms ┆ query_ms │ │ --- ┆ --- ┆ --- ┆ --- │ │ str ┆ str ┆ f64 ┆ f64 │ ╞════════════════════════╪═══════════════╪════════════╪══════════╡ │ GridInterpolator ┆ 1 point ┆ 255.055643 ┆ 0.072672 │ │ sax.interpolate_xarray ┆ 1 point ┆ 182.441207 ┆ 0.028726 │ │ GridInterpolator ┆ 10 000 points ┆ 289.342737 ┆ 0.336759 │ │ sax.interpolate_xarray ┆ 10 000 points ┆ 312.985474 ┆ 0.413262 │ └────────────────────────┴───────────────┴────────────┴──────────┘ ``` ## Symmetric CPW coupling: a three-dimensional dataset `qpdk/models/datasets/data/cpw_coupling.py` follows the same workflow, sweeping trace width, outer CPW slot width and the gap between two identical traces. The inner gap is fully etched; there is no ground strip between the traces. Both conductors and the outer ground rails define a two-dimensional transverse Palace mesh, eliminating end fringing. `Model.Lc` sets the implicit slice depth; stored capacitances are in F for that depth. The lookup has three geometry axes. ```bash uv run --script qpdk/models/datasets/data/cpw_coupling.py --dry-run uv run --script qpdk/models/datasets/data/cpw_coupling.py --sif /path/to/palace.sif uv run --script qpdk/models/datasets/data/cpw_coupling.py --topology as-drawn --sif /path/to/palace.sif ``` The generator refines every geometry until each raw Maxwell entry changes by at most 1% between successive meshes and identical traces agree within 1%. It stores the accepted finer matrix, mesh sizes, refinement level and measured change. This measures mesh sensitivity; it does not bound domain truncation or interpolation error. Separate domain, vacuum and held-out checks travel in the bundled dataset’s provenance. Solver matrices remain raw. The smaller development grid trades lookup resolution for faster regeneration; check its measured interpolation error before extending a design sweep. The following cells read the bundled Palace results and run no simulator. The same interpolator handles all three axes; fixing an axis only chooses a view of that three-dimensional lookup. ```ipython3 cpw_dataset = Dataset("cpw_coupling_palace") cpw_grid = cpw_dataset.grid("maxwell_capacitance") cpw_lookup = GridInterpolator(cpw_grid) slice_length = cpw_dataset.metadata.provenance["settings"]["slice_length_um"] * 1e-6 logger.info(f"CPW grid: {cpw_grid.values.shape}; domain: {cpw_lookup.domain}") ``` ```myst-ansi 2026-10-22 13:29:49.610 | INFO  | __main__::5 - CPW grid: (4, 4, 15, 2, 2); domain: {'width': (2.0, 30.0), 'cpw_gap': (1.0, 20.0), 'gap': (1.0, 250.0)} ``` ### Find simulations to refine Filter and project before collecting. The mesh change is the measured refinement check; the FEM norm ranks results to inspect and is not a capacitance error bound. Trace mismatch and comparison with conformal mapping provide independent flags. A flagged result needs investigation, even when its linear solve converged. Copy the CPW experiment beside the original, keep only flagged geometries in `GRID` or `GROUND_GRID`, and rerun with `--min-refinements 2`. This continues past an early small change. Use `--topology as-drawn` for the ground-strip grid. Keep the same `--workdir` and execution options to reuse matching completed solves. `--mesh-tolerance` and `--max-refinements` control the acceptance target and computation budget. ```ipython3 diagnostic_names = ( "fem_error_indicator_norm", "solver_relative_residual", "solver_iterations", "analytical_mutual_relative_difference", "mesh_relative_change", "mesh_refinement_level", "mesh_near_size", "mesh_far_size", ) diagnostics = ( cpw_dataset .scan() .filter(pl.col("quantity").is_in(diagnostic_names)) .select("width", "cpw_gap", "gap", "quantity", "value") .group_by("width", "cpw_gap", "gap") .agg( pl.col("value").filter(pl.col("quantity") == name).first().alias(name) for name in diagnostic_names ) ) symmetry = ( cpw_dataset .scan() .filter(pl.col("quantity") == "maxwell_capacitance", pl.col("row") == pl.col("col")) .select("width", "cpw_gap", "gap", "value") .group_by("width", "cpw_gap", "gap") .agg( trace_mismatch=(pl.col("value").max() - pl.col("value").min()) / pl.col("value").mean() ) ) indicator = pl.col("fem_error_indicator_norm") ( diagnostics .join(symmetry, on=["width", "cpw_gap", "gap"]) .with_columns( flag_for_refinement=( indicator > indicator.quantile(0.75) + 1.5 * (indicator.quantile(0.75) - indicator.quantile(0.25)) ) | (pl.col("mesh_relative_change") > 0.01) | (pl.col("trace_mismatch") > 0.01) | (pl.col("analytical_mutual_relative_difference").abs() > 0.03) ) .sort( pl.col("flag_for_refinement"), pl.col("analytical_mutual_relative_difference").abs(), pl.col("fem_error_indicator_norm"), descending=True, ) .head(10) .collect() ) ``` ```myst-ansi shape: (10, 13) ┌───────┬─────────┬───────┬─────────────┬───┬─────────────┬─────────────┬─────────────┬────────────┐ │ width ┆ cpw_gap ┆ gap ┆ fem_error_i ┆ … ┆ mesh_near_s ┆ mesh_far_si ┆ trace_misma ┆ flag_for_r │ │ --- ┆ --- ┆ --- ┆ ndicator_no ┆ ┆ ize ┆ ze ┆ tch ┆ efinement │ │ f64 ┆ f64 ┆ f64 ┆ rm ┆ ┆ --- ┆ --- ┆ --- ┆ --- │ │ ┆ ┆ ┆ --- ┆ ┆ f64 ┆ f64 ┆ f64 ┆ bool │ │ ┆ ┆ ┆ f64 ┆ ┆ ┆ ┆ ┆ │ ╞═══════╪═════════╪═══════╪═════════════╪═══╪═════════════╪═════════════╪═════════════╪════════════╡ │ 2.0 ┆ 1.0 ┆ 1.0 ┆ 0.202388 ┆ … ┆ 0.0000001 ┆ 0.00002 ┆ 0.000139 ┆ true │ │ 2.0 ┆ 6.0 ┆ 1.5 ┆ 0.188296 ┆ … ┆ 0.0000001 ┆ 0.00002 ┆ 0.000143 ┆ true │ │ 2.0 ┆ 1.0 ┆ 250.0 ┆ 0.142698 ┆ … ┆ 7.0000e-8 ┆ 0.000014 ┆ 0.000125 ┆ false │ │ 2.0 ┆ 20.0 ┆ 100.0 ┆ 0.144159 ┆ … ┆ 0.0000001 ┆ 0.00002 ┆ 0.001544 ┆ false │ │ 30.0 ┆ 1.0 ┆ 250.0 ┆ 0.119012 ┆ … ┆ 0.0000001 ┆ 0.00002 ┆ 0.000355 ┆ false │ │ 6.0 ┆ 1.0 ┆ 250.0 ┆ 0.138271 ┆ … ┆ 0.0000001 ┆ 0.00002 ┆ 0.000469 ┆ false │ │ 30.0 ┆ 3.0 ┆ 250.0 ┆ 0.11438 ┆ … ┆ 0.0000001 ┆ 0.00002 ┆ 0.00105 ┆ false │ │ 10.0 ┆ 1.0 ┆ 250.0 ┆ 0.107519 ┆ … ┆ 7.0000e-8 ┆ 0.000014 ┆ 0.000194 ┆ false │ │ 6.0 ┆ 1.0 ┆ 1.0 ┆ 0.166724 ┆ … ┆ 0.0000001 ┆ 0.00002 ┆ 0.000298 ┆ false │ │ 10.0 ┆ 3.0 ┆ 250.0 ┆ 0.12871 ┆ … ┆ 0.0000001 ┆ 0.00002 ┆ 0.000674 ┆ false │ └───────┴─────────┴───────┴─────────────┴───┴─────────────┴─────────────┴─────────────┴────────────┘ ``` ### A two-dimensional heatmap At a fixed outer slot width, vary trace width and inter-trace gap together. The black dots mark solved geometries; colours between them are interpolation. The logarithmic colour scale resolves weak coupling at large separations. ```ipython3 width_axis = jnp.linspace(*cpw_lookup.domain["width"], 101) gap_axis = jnp.geomspace(*cpw_lookup.domain["gap"], 151) W, G = jnp.meshgrid(width_axis, gap_axis, indexing="ij") mutual_per_length = -cpw_lookup(width=W, cpw_gap=6.0, gap=G)[..., 0, 1] / slice_length fig, ax = plt.subplots(constrained_layout=True) image = ax.pcolormesh( gap_axis, width_axis, mutual_per_length * 1e12, cmap=QPDK, norm=LogNorm(), shading="auto", rasterized=True, ) solved_width, solved_gap = jnp.meshgrid( cpw_grid.coords[0], cpw_grid.coords[2], indexing="ij" ) ax.scatter(solved_gap, solved_width, s=8, color="black") ax.set_xlabel(r"Inter-trace gap ($\text{µm}$)") ax.set_xscale("log") ax.set_ylabel(r"Trace width ($\text{µm}$)") ax.set_title(r"Outer slot width $6\,\text{µm}$") fig.colorbar(image, ax=ax, label=r"Mutual capacitance per length ($\text{pF/m}$)") plt.show() ``` ```myst-ansi
``` ### Slices through the three-dimensional lookup These three heatmaps vary all three parameters: the outer slot width changes between panels, while each panel scans trace width and inter-trace gap. A shared colour scale makes the effect of the third axis visible. The middle panel’s slot width lies between stored grid points. ```ipython3 outer_slots = jnp.array([3.0, 7.5, 12.0]) volume_slices = ( -cpw_lookup( width=W[None, ...], cpw_gap=outer_slots[:, None, None], gap=G[None, ...] )[..., 0, 1] / slice_length * 1e12 ) fig, axes = plt.subplots( 1, 3, figsize=(12, 3.8), sharex=True, sharey=True, constrained_layout=True ) heatmap_norm = LogNorm(vmin=float(volume_slices.min()), vmax=float(volume_slices.max())) for i, ax in enumerate(axes): image = ax.pcolormesh( gap_axis, width_axis, volume_slices[i], cmap=QPDK, norm=heatmap_norm, shading="auto", rasterized=True, ) ax.set_title(rf"Outer slot ${float(outer_slots[i]):g}\,\text{{µm}}$") ax.set_xlabel(r"Inter-trace gap ($\text{µm}$)") ax.set_xscale("log") axes[0].set_ylabel(r"Trace width ($\text{µm}$)") fig.colorbar( image, ax=list(axes), label=r"Mutual capacitance per length ($\text{pF/m}$)" ) plt.show() ``` ```myst-ansi
``` ### Compare with conformal mapping The analytical edge-coupled CPW expression assumes the same conductor sheets and dielectric half-spaces. Compare mutual capacitance per length directly, before turning either result into a circuit model. ```ipython3 analytical = cpw_cpw_coupling_capacitance_per_length_analytical( gap=gap_axis, width=10.0, cpw_gap=6.0, ep_r=cpw_dataset.metadata.provenance["settings"]["permittivity"], ) fem = -cpw_lookup(width=10.0, cpw_gap=6.0, gap=gap_axis)[..., 0, 1] / slice_length fig, axes = plt.subplots(1, 2, figsize=(10, 3.8), constrained_layout=True) axes[0].plot(gap_axis, fem * 1e12, label="Palace interpolation") axes[0].plot(gap_axis, analytical * 1e12, "--", label="Conformal mapping") axes[0].set_ylabel(r"Mutual capacitance per length ($\text{pF/m}$)") axes[0].legend() axes[1].plot(gap_axis, 100 * (fem / analytical - 1)) axes[1].set_ylabel(r"Difference from analytical ($\text{%}$)") for ax in axes: ax.set_xlabel(r"Inter-trace gap ($\text{µm}$)") ax.set_xscale("log") plt.show() ``` ```myst-ansi
``` ### Ground left between the CPW slots ![Fully etched inner gap compared with QPDK slots retaining a central ground strip.](notebooks/figures/cpw-coupling-ground.svg) `coupler_straight` etches a slot on each side of each trace. When the inter-trace gap exceeds twice the slot width, a ground strip remains between the slots. The `as-drawn` experiment keeps this strip at reference ground. It matches the fully etched geometry while the slots touch or overlap; beyond that point, the ground shields the traces. The FEM capacitances are solved independently; the analytical expression is used only for comparison. Both models share ideal sheets and dielectric half-spaces. The analytical ECCPW formula assumes a fully etched inner gap, so its difference from the shielded result measures a geometry change as well as numerical error. Positive strip widths are sampled logarithmically from 0.01 µm. Logarithmic interpolation follows the rapid shielding onset; narrower strips return NaN. ```ipython3 layout_dataset = Dataset("cpw_coupling_ground_strip_palace") layout_mutual = cpw_cpw_coupling_capacitance( f=5e9, length=1e6, gap=gap_axis, cross_section=coplanar_waveguide(width=10.0, gap=6.0), ) fig, ax = plt.subplots(figsize=(7, 4), constrained_layout=True) ax.loglog(gap_axis, fem * 1e12, label="Fully etched inner gap") ax.loglog(gap_axis, layout_mutual * 1e12, label="Drawn CPW slots") ax.axvline(12.0, color="gray", linestyle=":", label="Slots touch") ax.set_xlabel(r"Inter-trace gap ($\text{µm}$)") ax.set_ylabel(r"Mutual capacitance per length ($\text{pF/m}$)") ax.legend() plt.show() ``` ```myst-ansi
``` Compare both lookups with reference solves refined to a 0.1% successive-mesh change, larger domains, vacuum scaling and independent held-out geometries. Vacuum scaling uses byte-identical meshes to separate material scaling from refinement choices. These finite checks estimate sensitivity; they do not prove a uniform error bound. ```ipython3 quality_rows = [] for label, data in (("Etched gap", cpw_dataset), ("Ground strip", layout_dataset)): checks = data.metadata.provenance["validation"] quality_rows.append({ "Geometry": label, "Mesh (%)": 100 * checks["maximum_mesh_relative_change"], "Reference (%)": 100 * checks["representative_stricter_mesh_max_relative"], "Initial domain (%)": 100 * max( c["original_domain_relative"] for c in checks["mesh_domain"] if "original_domain_relative" in c ), "Domain (%)": 100 * max(c["domain_relative"] for c in checks["mesh_domain"]), "Same-mesh vacuum (%)": 100 * max(c["vacuum_scaling_relative"] for c in checks["mesh_domain"]), "Lookup (%)": 100 * checks["held_out_max_relative"], }) pl.DataFrame(quality_rows) ``` ```myst-ansi shape: (2, 7) ┌──────────────┬──────────┬───────────┬─────────────────────┬────────────┬────────────┬────────────┐ │ Geometry ┆ Mesh (%) ┆ Reference ┆ Initial domain (%) ┆ Domain (%) ┆ Same-mesh ┆ Lookup (%) │ │ --- ┆ --- ┆ (%) ┆ --- ┆ --- ┆ vacuum (%) ┆ --- │ │ str ┆ f64 ┆ --- ┆ f64 ┆ f64 ┆ --- ┆ f64 │ │ ┆ ┆ f64 ┆ ┆ ┆ f64 ┆ │ ╞══════════════╪══════════╪═══════════╪═════════════════════╪════════════╪════════════╪════════════╡ │ Etched gap ┆ 0.997168 ┆ 0.347286 ┆ 1.208178 ┆ 0.317497 ┆ 0.015447 ┆ 4.069779 │ │ Ground strip ┆ 0.99997 ┆ 0.78594 ┆ 1.27112 ┆ 0.173371 ┆ 0.003599 ┆ 4.478808 │ └──────────────┴──────────┴───────────┴─────────────────────┴────────────┴────────────┴────────────┘ ``` The widest gaps needed larger domains and finer meshes. Each correction records its raw solver source and the comparison with a separate finer reference. ```ipython3 pl.DataFrame([ { "Geometry": label, "Width (µm)": report["point"]["width"], "Slot (µm)": report["point"]["cpw_gap"], "Inner span (µm)": report["point"][ "gap" if label == "Etched gap" else "ground_strip_width" ], "Original pad (µm)": report["previous_domain_pad_um"], "Accepted pad (µm)": report["accepted_domain_pad_um"], "Original mesh": report["previous_refinement_level"], "Accepted mesh": report["accepted_refinement_level"], "Original/reference (%)": 100 * report["reference_max_relative_before"], "Accepted/reference (%)": 100 * report["reference_max_relative_after"], } for label, data in (("Etched gap", cpw_dataset), ("Ground strip", layout_dataset)) for report in data.metadata.provenance["curated_refinements"] ]) ``` ```myst-ansi shape: (2, 10) ┌────────────┬───────┬────────────┬────────────┬───┬───────────┬───────────┬───────────┬───────────┐ │ Geometry ┆ Width ┆ Slot (µm) ┆ Inner span ┆ … ┆ Original ┆ Accepted ┆ Original/ ┆ Accepted/ │ │ --- ┆ (µm) ┆ --- ┆ (µm) ┆ ┆ mesh ┆ mesh ┆ reference ┆ reference │ │ str ┆ --- ┆ f64 ┆ --- ┆ ┆ --- ┆ --- ┆ (%) ┆ (%) │ │ ┆ f64 ┆ ┆ f64 ┆ ┆ i64 ┆ i64 ┆ --- ┆ --- │ │ ┆ ┆ ┆ ┆ ┆ ┆ ┆ f64 ┆ f64 │ ╞════════════╪═══════╪════════════╪════════════╪═══╪═══════════╪═══════════╪═══════════╪═══════════╡ │ Etched gap ┆ 30.0 ┆ 20.0 ┆ 250.0 ┆ … ┆ 1 ┆ 5 ┆ 0.686527 ┆ 0.089547 │ │ Ground ┆ 30.0 ┆ 20.0 ┆ 248.0 ┆ … ┆ 1 ┆ 7 ┆ 0.84599 ┆ 0.068249 │ │ strip ┆ ┆ ┆ ┆ ┆ ┆ ┆ ┆ │ └────────────┴───────┴────────────┴────────────┴───┴───────────┴───────────┴───────────┴───────────┘ ``` ## A distributed four-port SAX model `cpw_cpw_coupling_capacitance()` uses the layout-matching lookup in `coupler_straight()`. For a distributed section, `cpw_coupling_model()` loads and validates the grid once, then returns a jittable model of the uniform coupled section. The slice capacitances determine even- and odd-mode impedances. Geometric inductance follows the quasi-TEM dielectric half-space relation; kinetic inductance, finite substrate thickness, dispersion and loss are omitted. Ports are `o1` lower-left, `o2` upper-left, `o3` upper-right and `o4` lower-right. Reference planes lie at the section ends. All four ports below use 50 ohm reference impedances. The remaining ports are matched when reading a single scattering entry; an open or short circuit must be connected explicitly in SAX. ```ipython3 cpw_model = cpw_coupling_model(layout_dataset) frequencies = jnp.linspace(0.1e9, 12e9, 301) fig, axes = plt.subplots(1, 2, figsize=(10, 3.8), constrained_layout=True) for gap in (2.0, 8.0, 25.0): scattering = jax.jit(cpw_model)( f=frequencies, length=1000.0, width=10.0, cpw_gap=6.0, gap=gap ) axes[0].plot( frequencies / 1e9, 20 * jnp.log10(jnp.abs(scattering["o1", "o2"])), label=rf"${gap:g}\,\text{{µm}}$", ) axes[1].plot( frequencies / 1e9, 20 * jnp.log10(jnp.abs(scattering["o1", "o4"])), label=rf"${gap:g}\,\text{{µm}}$", ) axes[0].set_ylabel(r"Near-end coupling $|S_{21}|$ ($\text{dB}$)") axes[1].set_ylabel(r"Through $|S_{41}|$ ($\text{dB}$)") for ax in axes: ax.set_xlabel(r"Frequency ($\text{GHz}$)") ax.legend(title="Inter-trace gap") plt.show() ``` ```myst-ansi
``` ## References