From 7f45d43b1ad1ac9be353a2edb1e00344b3e72c83 Mon Sep 17 00:00:00 2001 From: Yun Date: Mon, 7 Sep 2026 17:33:48 +0900 Subject: [PATCH 1/5] Evaluate the vacuum field once and reuse it across the time base The startup maps -- flux, poloidal field strength, decay index, the breakdown figure of merit -- are four readings of one field, so they are now four readings of one evaluation. Splitting them across separate computations is how they drift apart. The grid's response to each coil and vessel loop depends on geometry and not on time, so it is built once and cached; a different time then costs a contraction with that instant's currents. On the packaged shot that is 0.67 s once against 44 ms a frame, which is the difference between a time slider and a progress bar. Two details the cache had to get right. Its key comes from the same source extraction the matrices are built from, now factored out of compute_point_response_matrices_ods, so a cached matrix cannot outlive a geometry change it failed to notice. And each waveform is read once per frame rather than once per sample -- the smaller sibling of the defect fixed in #645. dpsi/dt is this map's own flux differenced against its neighbouring samples, so the induced electric field belongs to the same field as everything else instead of being spliced in from a measured loop voltage. The flux agrees with compute_null_ods to 0.0001% at the 95th percentile, which is what makes re-pointing the existing psi map safe. It disagrees only where a grid point lands on a source filament and both answers are meaningless, so the test bounds a percentile rather than a maximum. Note that compute_null_ods lays its grid out (Z, R) and this evaluator lays it out (R, Z). Co-Authored-By: Claude Opus 5 --- test/test_formula_docstrings.py | 1 + test/test_vacuum_field_map.py | 179 ++++++++++++++++++++++++ vaft/formula/equilibrium.py | 131 +++++++++++++++++ vaft/omas/process_wrapper.py | 239 +++++++++++++++++++++++++++++--- 4 files changed, 533 insertions(+), 17 deletions(-) create mode 100644 test/test_vacuum_field_map.py diff --git a/test/test_formula_docstrings.py b/test/test_formula_docstrings.py index df0341be..e151aeb0 100644 --- a/test/test_formula_docstrings.py +++ b/test/test_formula_docstrings.py @@ -45,6 +45,7 @@ "eK_from_K", "peaking_factor", "calculate_distance", + "poloidal_field_magnitude", "trapz_integral", "greens_integral_2d", "greens_integral_3d", diff --git a/test/test_vacuum_field_map.py b/test/test_vacuum_field_map.py new file mode 100644 index 00000000..140fb497 --- /dev/null +++ b/test/test_vacuum_field_map.py @@ -0,0 +1,179 @@ +"""The cached vacuum-field evaluator behind the interactive startup maps. + +One evaluation feeds every quantity the maps draw -- flux, |B_p|, the decay +index, the breakdown figure of merit -- so they cannot disagree with each other. +What has to hold is that the response matrices are built once per grid rather +than per frame, that the flux still matches the path it replaces, and that the +array orientation is the one the axes claim. +""" + +import numpy as np +import pytest + +import vaft.omas +from vaft.omas import process_wrapper as pw + + +@pytest.fixture(scope="module") +def solved(): + """The packaged shot with its vessel currents solved, once.""" + ods = vaft.omas.sample_ods() + vaft.omas.compute_eddy_currents(ods, [], []) + return ods + + +@pytest.fixture(autouse=True) +def _empty_cache(): + """No test may inherit another's cached matrices.""" + pw._VACUUM_MAP_CACHE.clear() + yield + pw._VACUUM_MAP_CACHE.clear() + + +COARSE = (np.linspace(0.12, 0.74, 9), np.linspace(-1.1, 1.1, 11)) + + +# --------------------------------------------------------------------------- +# Shape and orientation +# --------------------------------------------------------------------------- + +def test_every_field_is_indexed_r_then_z(solved): + """The axes are named; the arrays must actually be laid out that way. + + ``compute_null_ods`` returns the transpose of this, which is exactly the + kind of silent flip a plot inherits without complaint. + """ + result = pw.compute_vacuum_field_map(solved, time=0.29, grid=COARSE) + assert result["r"].size == 9 and result["z"].size == 11 + for key in ("psi", "b_r", "b_z", "dpsi_dt"): + assert result[key].shape == (9, 11), key + assert np.isfinite(result[key]).all(), key + + +def test_the_time_returned_is_a_stored_sample(solved): + """The map is snapped, not interpolated, so it must say where it landed.""" + time_base = np.asarray(solved["pf_active.time"], dtype=float) + result = pw.compute_vacuum_field_map(solved, time=0.2903177, grid=COARSE) + assert result["time"] == pytest.approx(time_base[result["time_index"]]) + assert abs(result["time"] - 0.2903177) <= np.diff(time_base).max() + + +# --------------------------------------------------------------------------- +# The cache, which is the whole reason a time slider is affordable +# --------------------------------------------------------------------------- + +def test_moving_only_in_time_does_not_rebuild_the_matrices(solved, monkeypatch): + """A slider move must cost a contraction, not a rebuild.""" + calls = [] + original = pw.compute_point_response_matrices_ods + + def counted(*args, **kwargs): + calls.append(1) + return original(*args, **kwargs) + + monkeypatch.setattr(pw, "compute_point_response_matrices_ods", counted) + first = pw.compute_vacuum_field_map(solved, time=0.29, grid=COARSE) + for t in (0.30, 0.31, 0.32): + pw.compute_vacuum_field_map(solved, time=t, grid=COARSE) + assert len(calls) == 1 + + # ... but a different grid is a different machine-to-grid response. + other = (COARSE[0], np.linspace(-1.0, 1.0, 11)) + pw.compute_vacuum_field_map(solved, time=0.29, grid=other) + assert len(calls) == 2 + + later = pw.compute_vacuum_field_map(solved, time=0.29, grid=COARSE) + assert len(calls) == 2 + np.testing.assert_allclose(later["psi"], first["psi"]) + + +def test_the_cache_is_keyed_on_the_geometry_it_describes(solved): + """Move a coil and the cached response no longer describes this machine.""" + before = pw._machine_fingerprint(solved) + original = float(solved["pf_active.coil.0.element.0.geometry.rectangle.r"]) + try: + solved["pf_active.coil.0.element.0.geometry.rectangle.r"] = original + 0.01 + assert pw._machine_fingerprint(solved) != before + finally: + solved["pf_active.coil.0.element.0.geometry.rectangle.r"] = original + assert pw._machine_fingerprint(solved) == before + + +def test_the_cache_does_not_grow_without_bound(solved): + """Each entry is three dense matrices; keeping every grid ever asked for + would be tens of gigabytes over a session.""" + for n in range(pw._VACUUM_MAP_CACHE_LIMIT + 2): + grid = (COARSE[0], np.linspace(-1.0 - 0.01 * n, 1.0, 7)) + pw.compute_vacuum_field_map(solved, time=0.29, grid=grid) + assert len(pw._VACUUM_MAP_CACHE) <= pw._VACUUM_MAP_CACHE_LIMIT + + +# --------------------------------------------------------------------------- +# Agreement with the path it replaces +# --------------------------------------------------------------------------- + +def test_the_flux_matches_compute_null_ods(solved): + """Re-pointing the vacuum psi map onto this evaluator must not move it. + + The tolerance is a percentile, not a maximum, and deliberately so: the two + paths differ only where a grid point sits on a source filament, where both + answers are meaningless. A max-based bound could not pass and would be + demanding the wrong thing. + """ + psi_reference, mesh_r, mesh_z = vaft.omas.compute_null_ods(solved, 0.29) + psi_reference = np.asarray(psi_reference, dtype=float) + r_axis, z_axis = np.unique(np.asarray(mesh_r)), np.unique(np.asarray(mesh_z)) + + result = pw.compute_vacuum_field_map(solved, time=0.29, grid=(r_axis, z_axis)) + # compute_null_ods lays its grid out (Z, R); this evaluator lays it out (R, Z). + reference = psi_reference.T + assert reference.shape == result["psi"].shape + + relative = np.abs(reference - result["psi"]) / np.maximum(np.abs(reference), 1e-12) + assert np.percentile(relative, 95) < 1e-4 + assert np.percentile(relative, 99) < 1e-3 + + +# --------------------------------------------------------------------------- +# The default grid and the magnitudes on it +# --------------------------------------------------------------------------- + +def test_the_default_grid_covers_the_limiter(solved): + """The map is read to site a null, so it covers where a plasma can sit.""" + wall_r = np.asarray(solved["wall.description_2d.0.limiter.unit.0.outline.r"], float) + wall_z = np.asarray(solved["wall.description_2d.0.limiter.unit.0.outline.z"], float) + r_axis, z_axis = pw._vacuum_map_grid(solved, 17) + assert r_axis.size == z_axis.size == 17 + assert r_axis[0] == pytest.approx(wall_r.min()) + assert r_axis[-1] == pytest.approx(wall_r.max()) + assert z_axis[0] == pytest.approx(wall_z.min()) + assert z_axis[-1] == pytest.approx(wall_z.max()) + + +def test_the_fields_are_the_size_a_vest_startup_is(solved): + """A guard on units and on the current assembly at once: a factor of 2*pi, + a sign flip on a coil or a wrong column order all leave this range.""" + result = pw.compute_vacuum_field_map(solved, time=0.29, grid=COARSE) + b_poloidal_gauss = np.hypot(result["b_r"], result["b_z"]) * 1e4 + assert 1.0 < np.median(b_poloidal_gauss) < 1000.0 + assert np.percentile(np.abs(result["psi"]), 50) < 1.0 # weber + assert np.any(np.abs(result["dpsi_dt"]) > 0.1) # a driven shot + + +def test_the_flux_derivative_is_the_slope_of_the_flux(solved): + """dpsi/dt has to be this map's own psi differenced, not a separate model.""" + time_base = np.asarray(solved["pf_active.time"], dtype=float) + middle = pw.compute_vacuum_field_map(solved, time=0.29, grid=COARSE) + index = middle["time_index"] + behind = pw.compute_vacuum_field_map(solved, time=time_base[index - 1], grid=COARSE) + ahead = pw.compute_vacuum_field_map(solved, time=time_base[index + 1], grid=COARSE) + expected = (ahead["psi"] - behind["psi"]) / (time_base[index + 1] - time_base[index - 1]) + np.testing.assert_allclose(middle["dpsi_dt"], expected, rtol=1e-10, atol=1e-12) + + +def test_a_geometry_only_sample_solves_its_own_eddy_currents(): + """compute_null_ods does this, so the evaluator replacing it must too.""" + ods = vaft.omas.sample_ods() + assert "time" not in ods["pf_passive"] + result = pw.compute_vacuum_field_map(ods, time=0.29, grid=COARSE) + assert np.isfinite(result["psi"]).all() diff --git a/vaft/formula/equilibrium.py b/vaft/formula/equilibrium.py index c6a96b5b..da4afa84 100644 --- a/vaft/formula/equilibrium.py +++ b/vaft/formula/equilibrium.py @@ -808,6 +808,137 @@ def bootstrap_current_fraction(n_e: float, # ------------------------------------------------------------------ +def poloidal_field_magnitude(b_r: np.ndarray, b_z: np.ndarray) -> np.ndarray: + r"""Poloidal field strength $|B_p| = \sqrt{B_R^2 + B_Z^2}$. + + Parameters + ---------- + b_r : array_like + Radial field component [T]. + b_z : array_like + Vertical field component [T]. + + Returns + ------- + np.ndarray + Poloidal field magnitude, elementwise [T]. + + Convention + ---------- + A magnitude, so it carries no COCOS sign: the orientation conventions cancel + in the quadrature. Both inputs must already be in tesla and in the same + convention as each other, which they are when they come from one call to + :func:`vaft.formula.green.green_br_bz_exact` or from one response matrix. + + Validity + -------- + Machine-independent. + """ + return np.hypot(np.asarray(b_r, dtype=float), np.asarray(b_z, dtype=float)) + + +def decay_index_from_bz( + r: np.ndarray, b_z: np.ndarray, *, axis: int = -1 +) -> np.ndarray: + r"""Field decay index $n = -\dfrac{R}{B_Z}\dfrac{\partial B_Z}{\partial R}$. + + How fast the vertical field falls off with major radius, which decides + whether the radial force balance holding a current ring is *stable*. + + Parameters + ---------- + r : array_like + Major radius, monotonic, along ``axis`` [m]. + b_z : array_like + Vertical field sampled on ``r``; may carry extra leading axes [T]. + axis : int, optional + Axis of ``b_z`` along which ``r`` varies [-]. + + Returns + ------- + np.ndarray + Decay index, ``nan`` where $B_Z$ vanishes [-]. + + Convention + ---------- + $0 < n < 1.5$ is the passively stable window: below zero the vertical field + does not restore a radial displacement, above 1.5 the ring is unstable to + vertical motion. + + Limitations + ----------- + Returns ``nan`` where $B_Z$ crosses zero rather than a large number. The + index is genuinely undefined on that surface, and a spike there is an + artefact of the division, not a physical instability -- which is what a + reader would otherwise take from it. + + Validity + -------- + Machine-independent. The stable window quoted above is the rigid-ring + result: it assumes a thin current ring and no conducting wall, so a real + vessel widens it. + + References + ---------- + .. [1] V. S. Mukhovatov and V. D. Shafranov, Nucl. Fusion 11 (1971) 605, + Sec. 2 (equilibrium of a current ring in a vertical field). + .. [2] J. Wesson, *Tokamaks*, 4th ed., Oxford University Press (2011), + Sec. 3.7 (vertical field and positional stability). + """ + r_arr = np.asarray(r, dtype=float) + b_arr = np.asarray(b_z, dtype=float) + gradient_bz = np.gradient(b_arr, r_arr, axis=axis, edge_order=2) + shape = [1] * b_arr.ndim + shape[axis] = r_arr.size + with np.errstate(divide="ignore", invalid="ignore"): + index = -r_arr.reshape(shape) * gradient_bz / b_arr + return np.where(np.isfinite(index), index, np.nan) + + +def toroidal_electric_field(r: np.ndarray, dpsi_dt: np.ndarray) -> np.ndarray: + r"""Toroidal electric field $E_\varphi = -\dfrac{1}{2\pi R}\dfrac{\partial\psi}{\partial t}$. + + The inductive drive a startup has to work with: the loop voltage + $-\partial\psi/\partial t$ spread around the torus at each major radius. + + Parameters + ---------- + r : array_like + Major radius [m]. + dpsi_dt : array_like + Time derivative of the **full-weber** poloidal flux, broadcastable + against ``r`` [Wb/s]. + + Returns + ------- + np.ndarray + Toroidal electric field [V/m]. + + Convention + ---------- + ``dpsi_dt`` is in weber, not weber per radian: the $2\pi$ here is the one + that turns a flux into a loop voltage, so passing a per-radian flux gives an + answer $2\pi$ too small. Green's-function flux + (:func:`vaft.formula.green.green_psi_exact`) is already full weber and needs + no conversion. + + Validity + -------- + Machine-independent. Faraday's law in axisymmetry, so it holds wherever the + flux does; it says nothing on its own about whether that field will break + the gas down, which also needs the connection length and the fill pressure. + + References + ---------- + .. [1] B. Lloyd et al., Nucl. Fusion 31 (1991) 2031, Sec. 2 (the toroidal + electric field required for tokamak start-up). + .. [2] J. Wesson, *Tokamaks*, 4th ed., Oxford University Press (2011), + Sec. 11.1 (start-up and breakdown). + """ + r_arr = np.asarray(r, dtype=float) + return -np.asarray(dpsi_dt, dtype=float) / (2.0 * np.pi * r_arr) + + def poloidal_field_factor( cocos: int | None, *, psi_per_radian: bool | None = None, ) -> float: diff --git a/vaft/omas/process_wrapper.py b/vaft/omas/process_wrapper.py index 110da2ef..6b90a105 100644 --- a/vaft/omas/process_wrapper.py +++ b/vaft/omas/process_wrapper.py @@ -326,6 +326,39 @@ def compute_grid_response_ods( logger.error(f"Error during computation: {e}") raise +def _vacuum_sources(ods: ODS) -> Tuple[List[float], List[float], List[float], List[int]]: + """Every current-carrying filament of the machine, as the field sees it. + + Coil elements enter individually with their signed turns; a passive loop + enters once, at its outline centroid (or its rectangle centre). ``groups`` + says which coil or loop each filament belongs to, so a response matrix can + be summed back down to one column per source. + + This is exactly what a vacuum response depends on besides the observation + points, which is why :func:`_machine_fingerprint` keys its cache on it. + """ + pf, pfp = ods["pf_active"], ods["pf_passive"] + nbcoil, nbloop = len(pf["coil"]), len(pfp["loop"]) + src_r, src_z, turns, groups = [], [], [], [] + for ii in range(nbcoil): + for jj in range(len(pf[f"coil.{ii}.element"])): + src_r.append(float(pf[f"coil.{ii}.element.{jj}.geometry.rectangle.r"])) + src_z.append(float(pf[f"coil.{ii}.element.{jj}.geometry.rectangle.z"])) + turns.append(float(pf[f"coil.{ii}.element.{jj}.turns_with_sign"])) + groups.append(ii) + for ii in range(nbloop): + geometry = pfp[f"loop.{ii}.element[0].geometry"] + if geometry["geometry_type"] == GEOMETRY_TYPE_POLYGON: + src_r.append(float(np.mean(geometry["outline.r"]))) + src_z.append(float(np.mean(geometry["outline.z"]))) + else: + src_r.append(float(geometry["rectangle.r"])) + src_z.append(float(geometry["rectangle.z"])) + turns.append(1.0) + groups.append(nbcoil + ii) + return src_r, src_z, turns, groups + + def compute_point_response_matrices_ods( ods: ODS, rz: List[List[float]], @@ -376,23 +409,7 @@ def compute_point_response_matrices_ods( nbcoil = len(pf["coil"]) nbloop = len(pfp["loop"]) - src_r, src_z, turns, groups = [], [], [], [] - for ii in range(nbcoil): - for jj in range(len(pf[f"coil.{ii}.element"])): - src_r.append(float(pf[f"coil.{ii}.element.{jj}.geometry.rectangle.r"])) - src_z.append(float(pf[f"coil.{ii}.element.{jj}.geometry.rectangle.z"])) - turns.append(float(pf[f"coil.{ii}.element.{jj}.turns_with_sign"])) - groups.append(ii) - for ii in range(nbloop): - geometry = pfp[f"loop.{ii}.element[0].geometry"] - if geometry["geometry_type"] == GEOMETRY_TYPE_POLYGON: - src_r.append(float(np.mean(geometry["outline.r"]))) - src_z.append(float(np.mean(geometry["outline.z"]))) - else: - src_r.append(float(geometry["rectangle.r"])) - src_z.append(float(geometry["rectangle.z"])) - turns.append(1.0) - groups.append(nbcoil + ii) + src_r, src_z, turns, groups = _vacuum_sources(ods) except KeyError as e: logger.error(f"Missing required data in ODS: {e}") raise @@ -859,6 +876,194 @@ def compute_point_vacuum_fields_ods( logger.error(f"Missing required data in ODS: {e}") raise +#: Response matrices for one grid and one machine, keyed by both. They are a +#: pure function of geometry, so a time slider pays for them once instead of on +#: every frame -- 2.4 s against 2.1 ms for the contraction that uses them. +_VACUUM_MAP_CACHE: dict[tuple, tuple] = {} + +#: How many cached grids to keep. Each is (n_points x n_sources) floats three +#: times over, ~97 MB at 65x65 on VEST, so this is deliberately small. +_VACUUM_MAP_CACHE_LIMIT = 4 + + +def _machine_fingerprint(ods: ODS) -> tuple: + """What the response matrices depend on, besides the observation points. + + Derived from :func:`_vacuum_sources` rather than restated, so a cached + matrix can never outlive a geometry change it did not notice. + """ + src_r, src_z, turns, groups = _vacuum_sources(ods) + return ( + np.asarray(src_r, dtype=float).tobytes(), + np.asarray(src_z, dtype=float).tobytes(), + np.asarray(turns, dtype=float).tobytes(), + tuple(groups), + ) + + +def _vacuum_response(ods: ODS, r_axis: ndarray, z_axis: ndarray): + """Cached ``(Psi, Bz, Br)`` response of one grid to every coil and loop.""" + key = ( + _machine_fingerprint(ods), + r_axis.tobytes(), + z_axis.tobytes(), + ) + cached = _VACUUM_MAP_CACHE.get(key) + if cached is None: + mesh_r, mesh_z = np.meshgrid(r_axis, z_axis, indexing="ij") + points = np.column_stack([mesh_r.ravel(), mesh_z.ravel()]).tolist() + cached = compute_point_response_matrices_ods(ods, points) + if len(_VACUUM_MAP_CACHE) >= _VACUUM_MAP_CACHE_LIMIT: + _VACUUM_MAP_CACHE.pop(next(iter(_VACUUM_MAP_CACHE))) + _VACUUM_MAP_CACHE[key] = cached + return cached + + +def _vacuum_currents(ods: ODS, indices) -> ndarray: + """Coil and passive-loop currents at the given samples of the PF time base. + + Returns ``(len(indices), nbcoil + nbloop)``, ordered ``[coils..., loops...]`` + to match the response matrices' columns. + + Every waveform is read once and then indexed, rather than re-read per + sample: on VEST that is 960 IMAS path parses instead of 960 per time, which + is the difference between a slider that moves and one that stutters. + """ + pf, pfp = ods["pf_active"], ods["pf_passive"] + take = np.asarray(indices, dtype=int).ravel() + columns = [np.asarray(pf[f"coil.{i}.current.data"], dtype=float)[take] + for i in range(len(pf["coil"]))] + columns += [np.asarray(pfp[f"loop.{i}.current"], dtype=float)[take] + for i in range(len(pfp["loop"]))] + return np.column_stack(columns) + + +def _vacuum_map_grid(ods: ODS, resolution: int) -> Tuple[ndarray, ndarray]: + """The default ``(R, Z)`` grid: the limiter's bounding box. + + The limiter is the region the field is being read for -- where breakdown + happens and where a null has to sit -- so that is what the map covers. It + is not source-free: on VEST several hundred coil filaments lie inside it, + and a grid point landing on one reads a meaningless spike. That is a + question for the contour levels, not for the extent. + """ + try: + wall_r = np.asarray( + ods["wall.description_2d.0.limiter.unit.0.outline.r"], dtype=float) + wall_z = np.asarray( + ods["wall.description_2d.0.limiter.unit.0.outline.z"], dtype=float) + except (KeyError, ValueError, IndexError): + wall_r = wall_z = np.asarray([]) + if wall_r.size < 3: + # No wall: fall back to the filaments' own extent, which at least + # brackets the machine. + src_r, src_z, _, _ = _vacuum_sources(ods) + wall_r, wall_z = np.asarray(src_r), np.asarray(src_z) + return ( + np.linspace(max(float(wall_r.min()), 1e-3), float(wall_r.max()), int(resolution)), + np.linspace(float(wall_z.min()), float(wall_z.max()), int(resolution)), + ) + + +def compute_vacuum_field_map( + ods: ODS, + *, + time: float, + grid: Optional[Tuple[ndarray, ndarray]] = None, + resolution: int = 65, +) -> Dict[str, ndarray]: + """Vacuum flux and poloidal field on an ``(R, Z)`` grid at one time. + + Everything the startup maps are drawn from -- the flux surfaces, |B_p|, the + decay index, the breakdown figure of merit -- comes from this one + evaluation, so they are guaranteed to describe the same field. + + The response of the grid to each coil and loop is a pure function of + geometry, so it is built once and cached; a different ``time`` then costs + only the contraction with that instant's currents. On VEST that is seconds + once against milliseconds a frame, which is what makes a time slider usable. + + Only vacuum sources contribute: coils and vessel eddy currents, no plasma. + Above roughly 10 kA of plasma current the map stops describing the machine. + + Parameters + ---------- + ods + Needs ``pf_active`` currents; passive-loop currents are solved with + :func:`compute_eddy_currents` if they are not already stored. + time + Seconds. Snapped to the nearest stored sample, not interpolated; the + sample actually used comes back in the result. + grid + ``(r_axis, z_axis)``. Defaults to the limiter's bounding box at + ``resolution`` points a side. + resolution + Points per axis for the default grid. Cost grows as its square: the + cached matrices are three arrays of ``resolution**2 x n_sources``. + + Returns + ------- + dict + ``r``, ``z`` (the axes), and ``psi`` [Wb], ``b_r``, ``b_z`` [T], + ``dpsi_dt`` [Wb/s], each shaped ``(len(r), len(z))``; plus ``time``, + the stored sample used, and ``time_index``, its position on the PF time + base. + + Notes + ----- + ``psi`` is full weber, as the Green's functions return it, not weber per + radian -- see :func:`vaft.data.eqdsk.ods_psi_to_wb_per_radian_factor`. + """ + if "time" not in ods["pf_passive"]: + # Geometry-only samples carry no passive-loop waveform; solve it rather + # than refusing, matching compute_null_ods. + compute_eddy_currents(ods, plasma=[], ip=[]) + + if grid is None: + r_axis, z_axis = _vacuum_map_grid(ods, resolution) + else: + r_axis = np.asarray(grid[0], dtype=float).ravel() + z_axis = np.asarray(grid[1], dtype=float).ravel() + if r_axis.size == 0 or z_axis.size == 0: + raise ValueError("grid axes must each hold at least one point") + + time_base = np.asarray(ods["pf_active.time"], dtype=float) + index = int(np.argmin(np.abs(time_base - float(time)))) + + psi_response, bz_response, br_response = _vacuum_response(ods, r_axis, z_axis) + width = len(ods["pf_active.coil"]) + len(ods["pf_passive.loop"]) + psi_response = psi_response[:, :width] + shape = (r_axis.size, z_axis.size) + + # dpsi/dt comes from the model's own flux at the neighbouring samples, so + # the induced electric field belongs to the same field as everything else + # rather than being spliced in from a measured loop voltage. All three + # samples are read in one pass over the waveforms. + lo, hi = max(index - 1, 0), min(index + 1, time_base.size - 1) + behind, currents, ahead = _vacuum_currents(ods, (lo, index, hi)) + + psi = (psi_response @ currents).reshape(shape) + b_z = (bz_response[:, :width] @ currents).reshape(shape) + b_r = (br_response[:, :width] @ currents).reshape(shape) + + span = float(time_base[hi] - time_base[lo]) + if span > 0.0: + dpsi_dt = ((psi_response @ (ahead - behind)) / span).reshape(shape) + else: + dpsi_dt = np.zeros(shape) + + return { + "r": r_axis, + "z": z_axis, + "psi": psi, + "b_r": b_r, + "b_z": b_z, + "dpsi_dt": dpsi_dt, + "time": float(time_base[index]), + "time_index": index, + } + + def compute_null_ods(ods, time): """Compute poloidal flux (psi) on grid at given time using coil and eddy currents. From 28096845edc729b012dec55d4be4941c1bcf7d11 Mon Sep 17 00:00:00 2001 From: Yun Date: Mon, 7 Sep 2026 17:54:47 +0900 Subject: [PATCH 2/5] Draw the vacuum field a start-up is actually read from Four quantities on one plot -- the flux, |B_p| in gauss, the field decay index, and the breakdown figure E_t B_t / B_p -- chosen with field=, and stepped along the PF time base with time_index=. They come from one cached evaluation, so the panels of a start-up figure cannot end up describing four slightly different instants. The map is confined to the limiter outline, and that is the load-bearing choice. VEST's limiter bounding box contains 330 of the machine's 1480 current filaments, so a map over the box is dominated by conductors: at the breakdown instant the whole vessel fell in the bottom contour bin while a handful of coil points took the entire colour range. Not one filament lies inside the limiter outline itself, so masking to it leaves exactly the field a discharge would see -- and the levels then describe that field. Levels come from a percentile of what is left, and Field2D gained an extend policy so the excluded points saturate at the end colours instead of coming out as holes indistinguishable from no data. A dense time base is a slider, not a list. This is the first plot to emit a range control, declared through a times block that names the option it drives, so a stored equilibrium's handful of slices keeps its radio buttons. Discovery resolves the default instant the way the builder does -- from the magnetics, before the private copy that drops them -- so opening the controls does not move the figure. Units are three vocabularies here, not one: gauss for the field, V/m for the breakdown figure, nothing at all for the decay index. The existing applies_to gate handles it, offering the flux units only while the flux is drawn. The vessel currents are cached alongside the response matrices. A recipe hands the evaluator a private copy each frame, which would otherwise repeat the same two-second eddy solve every time; a frame now costs 0.22 s from a caller's untouched ODS. The plan also had equilibrium_field_psi_vacuum move onto this evaluator. Measurement says no: on the 129x129 equilibrium grid that plot draws, the exact elliptic kernels cost 18 s against the polynomial path's 5.8 s and their cached matrices would be 591 MB. That trade buys a slider, not a single figure. It stays where it is, and the evaluator now refuses a grid whose response would not fit its budget rather than quietly allocating half a gigabyte and keeping it. Co-Authored-By: Claude Opus 5 --- test/test_vacuum_field_plot.py | 180 +++++++++++++++++++++++++ vaft/imas/plotting.py | 19 +++ vaft/omas/plotting.py | 21 +++ vaft/omas/process_wrapper.py | 73 ++++++++++- vaft/plot/__init__.py | 1 + vaft/plot/backend/discovery.py | 66 +++++++++- vaft/plot/backend/options.py | 5 + vaft/plot/backend/recipes.py | 231 +++++++++++++++++++++++++++++++++ vaft/plot/controls.py | 21 ++- vaft/plot/display.py | 10 +- vaft/plot/models.py | 11 ++ vaft/plot/renderers/fields.py | 24 ++++ vaft/plot/taxonomy.py | 3 + 13 files changed, 659 insertions(+), 6 deletions(-) create mode 100644 test/test_vacuum_field_plot.py diff --git a/test/test_vacuum_field_plot.py b/test/test_vacuum_field_plot.py new file mode 100644 index 00000000..14bbec96 --- /dev/null +++ b/test/test_vacuum_field_plot.py @@ -0,0 +1,180 @@ +"""The interactive vacuum-field map: what it draws and what it offers to change. + +Four quantities -- flux, |B_p|, the decay index, the breakdown figure of merit +-- from one cached evaluation, over the PF time base rather than a handful of +stored equilibrium slices. +""" + +import numpy as np +import pytest + +import vaft.omas +from vaft.omas import process_wrapper as pw +from vaft.plot.backend import recipes +from vaft.plot.backend.discovery import describe_one +from vaft.plot.controls import controls_for +from vaft.plot.models import Field2D + +COARSE = 21 + + +@pytest.fixture(scope="module") +def ods(): + return vaft.omas.sample_ods() + + +@pytest.fixture(scope="module") +def record(ods): + return describe_one("vacuum_field", [("sample", ods)]) + + +# --------------------------------------------------------------------------- +# What it draws +# --------------------------------------------------------------------------- + +@pytest.mark.parametrize("field", recipes.VACUUM_FIELD_NAMES) +def test_each_quantity_builds_a_field_on_the_poloidal_plane(ods, field): + model = recipes._build_vacuum_field(ods, field=field, resolution=COARSE) + assert isinstance(model, Field2D) + assert model.values.shape == (model.z.size, model.r.size) + assert np.isfinite(model.values).any() + + +def test_the_quantities_are_read_off_one_field(ods): + """|B_p| must be the quadrature of the same B_r and B_z the map holds -- + four separate computations is how the panels of a startup figure end up + describing four slightly different instants.""" + model = recipes._build_vacuum_field(ods, field="b_poloidal", resolution=COARSE) + instant = recipes._vacuum_psi_time(ods) + result = pw.compute_vacuum_field_map( + recipes._isolated_copy(ods, recipes._NULL_FIELD_ROOTS), + time=instant, + grid=(model.r, model.z), + ) + expected = np.hypot(result["b_r"], result["b_z"]).T * 1e4 # gauss + inside = np.isfinite(model.values) + np.testing.assert_allclose(model.values[inside], expected[inside], rtol=1e-9) + + +def test_a_misspelt_quantity_names_the_vocabulary(ods): + with pytest.raises(ValueError, match="b_poloidal"): + recipes._build_vacuum_field(ods, field="b_pol") + + +def test_the_map_is_confined_to_the_limiter(ods): + """Every current filament of the machine lies outside the limiter outline, + so the plasma-facing region is exactly the part of a vacuum map that is not + a conductor's own singular field.""" + model = recipes._build_vacuum_field(ods, field="b_poloidal", resolution=COARSE) + drawn = np.isfinite(model.values) + assert drawn.any() and not drawn.all() + + source_r, source_z, _, _ = pw._vacuum_sources(ods) + interior = recipes._limiter_interior(ods, model.r, model.z) + assert interior is not None + from matplotlib.path import Path + + outline = recipes._wall_layers(ods)[0] + polygon = Path(np.column_stack([outline.r, outline.z])) + assert not polygon.contains_points(np.column_stack([source_r, source_z])).any() + + +def test_the_units_are_the_ones_a_start_up_is_read_in(ods): + """Gauss for the poloidal field, V/m for the breakdown figure, nothing at + all for the decay index -- three vocabularies, one per quantity.""" + units = { + field: recipes._build_vacuum_field(ods, field=field, resolution=COARSE).display.unit + for field in recipes.VACUUM_FIELD_NAMES + } + assert units == {"psi": "mWb", "b_poloidal": "G", "decay_index": "", "breakdown": "V/m"} + + +def test_out_of_range_values_saturate_rather_than_vanish(ods): + """The levels come from a percentile, so the points it excludes must come + out as the end colour and not as holes indistinguishable from no data.""" + for field in ("b_poloidal", "breakdown", "decay_index"): + model = recipes._build_vacuum_field(ods, field=field, resolution=COARSE) + assert model.extend in ("max", "both"), field + assert recipes._build_vacuum_field(ods, field="psi", resolution=COARSE).extend == "neither" + + +def test_the_decay_index_marks_its_stable_band(ods): + model = recipes._build_vacuum_field(ods, field="decay_index", resolution=COARSE) + assert model.secondary_levels == recipes.DECAY_INDEX_STABLE_BAND + + +# --------------------------------------------------------------------------- +# Which instant +# --------------------------------------------------------------------------- + +def test_time_index_selects_a_pf_sample(ods): + base = np.asarray(ods["pf_active.time"], dtype=float) + model = recipes._build_vacuum_field(ods, field="psi", time_index=1200, resolution=COARSE) + assert f"{base[1200] * 1e3:.1f} ms" in model.title + + +def test_an_out_of_range_time_index_says_how_many_samples_there_are(ods): + with pytest.raises(ValueError, match="stored PF samples"): + recipes._build_vacuum_field(ods, field="psi", time_index=99999, resolution=COARSE) + + +def test_the_slider_starts_where_the_plot_does(ods, record): + """The default instant is the breakdown onset, read from the magnetics. + Discovery has to resolve it the same way the builder does, or opening the + controls jumps the figure to a different time.""" + default_index = record.times["selected"] + base = np.asarray(ods["pf_active.time"], dtype=float) + drawn = recipes._build_vacuum_field(ods, field="psi", resolution=COARSE) + assert f"{base[default_index] * 1e3:.1f} ms" in drawn.title + + +# --------------------------------------------------------------------------- +# What it offers to change +# --------------------------------------------------------------------------- + +def test_a_dense_time_base_is_offered_as_a_slider(record): + """Thousands of PF samples cannot be radio buttons, which is what a stored + equilibrium's handful of slices gets.""" + controls = {control.name: control for control in controls_for(record)} + assert "time_slice" not in controls + slider = controls["time_index"] + assert slider.kind == "range" + assert slider.options == (0, record.times["count"] - 1, 1) + assert slider.options[0] <= slider.default <= slider.options[1] + assert "ms" in slider.label + + +def test_the_quantity_is_offered_with_its_own_vocabulary(record): + controls = {control.name: control for control in controls_for(record)} + assert controls["field"].options == recipes.VACUUM_FIELD_NAMES + assert controls["field"].default == "psi" + + +def test_the_flux_units_are_offered_only_while_the_flux_is_drawn(record): + """Gauss, V/m and dimensionless are three different vocabularies; a unit + control that ignored which quantity is drawn would offer mWb for |B_p|.""" + controls = {control.name: control for control in controls_for(record)} + assert controls["units"].applies_to == {"field": ("psi",)} + + +def test_a_vacuum_map_offers_no_flux_map_style(record): + """The psi styles normalise against a separatrix this map does not have.""" + assert "style" not in {control.name for control in controls_for(record)} + + +def test_the_breakdown_figure_is_withheld_without_a_toroidal_field(ods): + """A control that would raise when used is worse than one that is absent.""" + stripped = recipes._isolated_copy(ods, ("pf_active", "pf_passive", "wall", "equilibrium")) + record = describe_one("vacuum_field", [("no tf", stripped)]) + assert "breakdown" not in record.fields["options"] + assert set(record.fields["options"]) == {"psi", "b_poloidal", "decay_index"} + with pytest.raises(ValueError, match="tf.b_field_tor_vacuum_r"): + recipes._build_vacuum_field(stripped, field="breakdown", resolution=COARSE) + + +def test_the_caller_s_ods_is_left_as_it_was_found(ods): + """The evaluator solves the vessel currents; that must not land in the + caller's ODS.""" + assert "time" not in ods["pf_passive"] + recipes._build_vacuum_field(ods, field="psi", resolution=COARSE) + assert "time" not in ods["pf_passive"] diff --git a/vaft/imas/plotting.py b/vaft/imas/plotting.py index 1e1c6168..f297ab25 100644 --- a/vaft/imas/plotting.py +++ b/vaft/imas/plotting.py @@ -571,6 +571,24 @@ def plot_equilibrium_field_psi_vacuum( return render("equilibrium_field_psi_vacuum", source, ax=ax, show=show, label=label, **options) +def plot_vacuum_field( + source: Any, + *, + ax: Any = None, + show: bool = False, + label: str | Sequence[str] = "shot", + **options: Any, +) -> tuple[Any, Any]: + """One quantity of the coils' and vessel's vacuum field, at one instant. + + ``field=`` chooses among ``psi``, ``b_poloidal``, ``decay_index`` and + ``breakdown``; ``time_index=`` steps along the PF time base. + + Renders with :func:`vaft.plot.vacuum_field` from native IMAS input. + """ + return render("vacuum_field", source, ax=ax, show=show, label=label, **options) + + def plot_equilibrium_geometry_boundary( source: Any, *, @@ -1997,6 +2015,7 @@ def plot_wall_geometry_poloidal( "plot_equilibrium_field_2d", "plot_equilibrium_field_psi", "plot_equilibrium_field_psi_vacuum", + "plot_vacuum_field", "plot_equilibrium_geometry_boundary", "plot_equilibrium_geometry_topview", "plot_equilibrium_overview", diff --git a/vaft/omas/plotting.py b/vaft/omas/plotting.py index d89afecc..92aadd86 100644 --- a/vaft/omas/plotting.py +++ b/vaft/omas/plotting.py @@ -807,6 +807,26 @@ def plot_equilibrium_field_psi_vacuum( ) +def plot_vacuum_field( + source: Any, + *, + ax: Any = None, + show: bool = False, + label: str | Sequence[str] = "shot", + **options: Any, +) -> tuple[Any, Any]: + """One quantity of the coils' and vessel's vacuum field, at one instant. + + ``field=`` chooses among ``psi``, ``b_poloidal``, ``decay_index`` and + ``breakdown``; ``time_index=`` steps along the PF time base. + + Renders with :func:`vaft.plot.vacuum_field`. + """ + return render( + "vacuum_field", source, ax=ax, show=show, label=label, **options + ) + + def plot_equilibrium_geometry_boundary( source: Any, *, @@ -2777,6 +2797,7 @@ def plot_nbi_profile_current_drive( "plot_equilibrium_field_2d", "plot_equilibrium_field_psi", "plot_equilibrium_field_psi_vacuum", + "plot_vacuum_field", "plot_equilibrium_geometry_boundary", "plot_equilibrium_geometry_topview", "plot_equilibrium_overview", diff --git a/vaft/omas/process_wrapper.py b/vaft/omas/process_wrapper.py index 6b90a105..bb9672d3 100644 --- a/vaft/omas/process_wrapper.py +++ b/vaft/omas/process_wrapper.py @@ -885,6 +885,15 @@ def compute_point_vacuum_fields_ods( #: times over, ~97 MB at 65x65 on VEST, so this is deliberately small. _VACUUM_MAP_CACHE_LIMIT = 4 +#: Vessel currents solved from one PF programme, keyed by that programme and +#: the machine. A plot that hands the evaluator a private copy of the ODS each +#: time -- which is how a recipe avoids mutating its caller's -- would otherwise +#: repeat the same two-second eddy solve on every frame. +_VACUUM_EDDY_CACHE: dict[tuple, tuple] = {} + +#: Two entries: each holds one current per loop per sample, ~19 MB on VEST. +_VACUUM_EDDY_CACHE_LIMIT = 2 + def _machine_fingerprint(ods: ODS) -> tuple: """What the response matrices depend on, besides the observation points. @@ -919,6 +928,42 @@ def _vacuum_response(ods: ODS, r_axis: ndarray, z_axis: ndarray): return cached +def _solve_or_recall_eddy_currents(ods: ODS) -> None: + """Fill in ``pf_passive`` currents, reusing an identical earlier solve. + + The vessel's response is fixed by the machine and the PF programme, so two + ODS objects carrying the same coil waveforms have the same answer. Solving + it again costs about two seconds, which is most of what an interactive + frame would otherwise spend. + """ + pf = ods["pf_active"] + time_base = np.asarray(pf["time"], dtype=float) + coil_currents = np.column_stack([ + np.asarray(pf[f"coil.{i}.current.data"], dtype=float) + for i in range(len(pf["coil"])) + ]) if len(pf["coil"]) else np.zeros((time_base.size, 0)) + key = (_machine_fingerprint(ods), time_base.tobytes(), coil_currents.tobytes()) + + cached = _VACUUM_EDDY_CACHE.get(key) + if cached is None: + compute_eddy_currents(ods, plasma=[], ip=[]) + pfp = ods["pf_passive"] + cached = ( + np.asarray(pfp["time"], dtype=float), + np.column_stack([np.asarray(pfp[f"loop.{i}.current"], dtype=float) + for i in range(len(pfp["loop"]))]), + ) + if len(_VACUUM_EDDY_CACHE) >= _VACUUM_EDDY_CACHE_LIMIT: + _VACUUM_EDDY_CACHE.pop(next(iter(_VACUUM_EDDY_CACHE))) + _VACUUM_EDDY_CACHE[key] = cached + return + + passive_time, loop_currents = cached + ods["pf_passive.time"] = passive_time + for index in range(loop_currents.shape[1]): + ods[f"pf_passive.loop.{index}.current"] = loop_currents[:, index] + + def _vacuum_currents(ods: ODS, indices) -> ndarray: """Coil and passive-loop currents at the given samples of the PF time base. @@ -938,6 +983,31 @@ def _vacuum_currents(ods: ODS, indices) -> ndarray: return np.column_stack(columns) +#: How much one cached grid's response may occupy, in bytes. Three dense +#: (n_points x n_sources) float64 matrices, so it grows with the fourth power +#: of the resolution -- 97 MB at 65 on VEST, 591 MB at 129. +_VACUUM_RESPONSE_BUDGET = 200 * 1024 ** 2 + + +def _refuse_oversized_grid(ods: ODS, r_axis: ndarray, z_axis: ndarray) -> None: + """Refuse a grid whose cached response would not fit, naming one that does. + + Silently allocating several hundred megabytes -- and then keeping it, since + the point of the cache is to keep it -- is worse than saying no. + """ + sources = len(_vacuum_sources(ods)[0]) + needed = 3 * r_axis.size * z_axis.size * sources * 8 + if needed <= _VACUUM_RESPONSE_BUDGET: + return + affordable = int(np.sqrt(_VACUUM_RESPONSE_BUDGET / (3 * sources * 8))) + raise ValueError( + f"a {r_axis.size}x{z_axis.size} grid over {sources} source filaments " + f"needs {needed / 1024 ** 2:.0f} MB of cached response matrices, above " + f"the {_VACUUM_RESPONSE_BUDGET / 1024 ** 2:.0f} MB budget; " + f"resolution={affordable} fits" + ) + + def _vacuum_map_grid(ods: ODS, resolution: int) -> Tuple[ndarray, ndarray]: """The default ``(R, Z)`` grid: the limiter's bounding box. @@ -1017,7 +1087,7 @@ def compute_vacuum_field_map( if "time" not in ods["pf_passive"]: # Geometry-only samples carry no passive-loop waveform; solve it rather # than refusing, matching compute_null_ods. - compute_eddy_currents(ods, plasma=[], ip=[]) + _solve_or_recall_eddy_currents(ods) if grid is None: r_axis, z_axis = _vacuum_map_grid(ods, resolution) @@ -1030,6 +1100,7 @@ def compute_vacuum_field_map( time_base = np.asarray(ods["pf_active.time"], dtype=float) index = int(np.argmin(np.abs(time_base - float(time)))) + _refuse_oversized_grid(ods, r_axis, z_axis) psi_response, bz_response, br_response = _vacuum_response(ods, r_axis, z_axis) width = len(ods["pf_active.coil"]) + len(ods["pf_passive.loop"]) psi_response = psi_response[:, :width] diff --git a/vaft/plot/__init__.py b/vaft/plot/__init__.py index 6f463689..85a2ef35 100755 --- a/vaft/plot/__init__.py +++ b/vaft/plot/__init__.py @@ -204,6 +204,7 @@ def plasma_current_time(model: LineSeries, *, ax=None, show=False, **style equilibrium_field_psi, equilibrium_field_psi_vacuum, passive_structure_field_wall_reduction, + vacuum_field, ) from .renderers.geometry import ( charge_exchange_geometry_poloidal, diff --git a/vaft/plot/backend/discovery.py b/vaft/plot/backend/discovery.py index a3ca08ed..7c2a4e4d 100644 --- a/vaft/plot/backend/discovery.py +++ b/vaft/plot/backend/discovery.py @@ -286,14 +286,19 @@ def _declare(record: PlotCapability) -> PlotCapability: "options": tuple(recipe.coordinates), "declared": tuple(recipe.coordinates), } - elif record.name in PSI_FIELD_CONVENTIONS: + elif record.name in PSI_FIELD_CONVENTIONS or record.name == "vacuum_field": # The psi maps take units= (issue #478); the default unit follows the # stored convention, which only an input can tell -- see _evaluate. + # The vacuum map is computed from the Green's functions, so its flux + # is full weber whatever the stored equilibrium declares, and it can + # say so without an input. + convention = PSI_FIELD_CONVENTIONS.get(record.name, "Wb") updates["display"] = { - "unit": None, + "unit": resolve_display(convention, subject="equilibrium").unit + if convention else None, "units": allowed_units("magnetic_flux", "equilibrium"), "notation": "auto", - "convention": PSI_FIELD_CONVENTIONS[record.name], + "convention": convention, } if isinstance(recipe, LineRecipe): options = abscissa_options(recipe) @@ -428,6 +433,8 @@ def _evaluate(record: PlotCapability, entries: Sequence[tuple[str, Any]]) -> Plo updates["coordinates"] = _coordinates_block(record, recipe, ods) if _takes_time_slice(record.name): updates["slices"] = _slices_block(ods) + if record.name == "vacuum_field": + updates["times"] = _pf_samples_block(ods) if record.fields: updates["fields"] = _fields_block(record, ods) if record.name in PSI_FIELD_CONVENTIONS: @@ -472,6 +479,8 @@ def _fields_block(record: PlotCapability, ods: Any) -> dict[str, Any]: from .recipes import _count, resolve_time_slice declared = tuple(record.fields.get("declared") or record.fields.get("options") or ()) + if record.name == "vacuum_field": + return _vacuum_fields_block(record, ods, declared) if not _count(ods, "equilibrium.time_slice"): return {**record.fields, "options": ()} try: @@ -492,6 +501,28 @@ def _fields_block(record: PlotCapability, ods: Any) -> dict[str, Any]: } +def _vacuum_fields_block( + record: PlotCapability, ods: Any, declared: tuple[str, ...] +) -> dict[str, Any]: + """The vacuum quantities this input can draw. + + The field itself needs only the PF programme, which the plot already + requires; the breakdown figure of merit additionally needs the toroidal + field, so an input without it is offered the other three rather than a + control that would raise. + """ + options = tuple( + name for name in declared + if name != "breakdown" or _has(ods, "tf.b_field_tor_vacuum_r.data") + ) + default = record.fields.get("default") + return { + "default": default if default in options else (options[0] if options else None), + "options": options, + "declared": declared, + } + + def _coordinates_block(record: PlotCapability, recipe: ProfileRecipe, ods: Any) -> dict[str, Any]: """The coordinates this input can resolve for a slice-indexed profile. @@ -705,6 +736,35 @@ def _times_block(ods: Any, recipe: Any) -> dict[str, Any]: return {} +def _pf_samples_block(ods: Any) -> dict[str, Any]: + """The PF time base a vacuum map steps along, as a dense indexed axis. + + ``option`` names the keyword the index is passed as, which is what tells + the control layer to offer a slider rather than a list: an equilibrium has + a handful of stored slices to pick between, a PF programme has thousands of + samples. ``selected`` is the sample the plot draws on its own, so opening + the controls does not move the figure. + """ + from .recipes import _array, _vacuum_map_time + + axis = _array(ods, "pf_active.time") + if axis is None or len(axis) < 2: + return {} + axis = np.asarray(axis, dtype=float) + try: + default = _vacuum_map_time(ods, None, None) + selected = int(np.argmin(np.abs(axis - float(default)))) + except (ValueError, KeyError, IndexError): + selected = 0 + return { + "start": float(axis[0]), + "stop": float(axis[-1]), + "count": int(axis.size), + "option": "time_index", + "selected": selected, + } + + def _validity_block(*, present: bool, flagged: int) -> dict[str, Any]: if not present: return {} diff --git a/vaft/plot/backend/options.py b/vaft/plot/backend/options.py index 8ce23662..d4182c09 100644 --- a/vaft/plot/backend/options.py +++ b/vaft/plot/backend/options.py @@ -67,6 +67,11 @@ def _specs() -> tuple[OptionSpec, ...]: OptionSpec("time_slice", "int", description="stored equilibrium slice index"), OptionSpec("time", "float", description="a time in seconds, snapped to a stored slice"), OptionSpec("time_range", "range", description="(start, stop) in seconds"), + # A dense time base is indexed, not chosen from a list: the vacuum map + # runs over the PF samples, thousands of them, where time_slice= names + # one of a handful of stored equilibria. + OptionSpec("time_index", "int", description="position on a dense time base"), + OptionSpec("resolution", "int", description="points per axis of a computed 2-D grid"), OptionSpec("centre", "range", description="(r0, z0) in metres the poloidal angle is measured about"), OptionSpec("angle", "choice", "recipes.ANGLE_SOURCES", "where a sensor's poloidal angle comes from"), OptionSpec("overlay", "multi", "recipes.CAMERA_OVERLAYS", "what is drawn over a map"), diff --git a/vaft/plot/backend/recipes.py b/vaft/plot/backend/recipes.py index 7228704f..f32dec8d 100644 --- a/vaft/plot/backend/recipes.py +++ b/vaft/plot/backend/recipes.py @@ -3009,6 +3009,14 @@ def _build_vacuum_psi( ) -> Field2D: """Vacuum poloidal flux from the PF coils, via the OMAS null-field helper. + Stays on the polynomial-kernel path rather than the cached evaluator behind + ``vacuum_field``. Measured on the packaged shot: this plot draws on the + equilibrium's own 129x129 grid, where the exact elliptic Green's functions + cost 18 s against 5.8 s and their cached response matrices would be 591 MB. + That trade pays for a slider and not for a single figure, which is what + this plot is; the two agree on psi to 0.0002% at the 95th percentile + either way. + The Green's functions give full weber, whatever convention the stored equilibrium uses, so the map is labelled from ``"Wb"`` through the display policy (mWb by default; ``units=`` chooses, per-radian included). @@ -3040,6 +3048,220 @@ def _build_vacuum_psi( ) +@dataclass(frozen=True) +class VacuumField: + """One quantity the vacuum map can draw, and how it is displayed. + + ``canonical_unit`` is empty for the decay index, which reaches the display + policy through :data:`vaft.plot.display.DIMENSIONLESS_DISPLAY` instead. + ``upper`` is the percentile the default contour levels stop at: the field + diverges on the coil filaments the grid sits among, so raw min-to-max + levels would spend the whole colour range on a handful of points and leave + the vessel blank. + """ + + name: str + label: str + canonical_unit: str + subject: str = "vacuum" + filled: bool = True + upper: float = 95.0 + lower: float | None = None + levels: int = 24 + + @property + def extend(self) -> str: + """How the map saturates outside its levels -- see :class:`Field2D`.""" + if not self.filled: + return "neither" + return "max" if self.lower is not None else "both" + + +VACUUM_FIELDS: dict[str, VacuumField] = { + field.name: field + for field in ( + VacuumField("psi", "Vacuum Poloidal Flux", "Wb", subject="equilibrium", + filled=False, upper=100.0, levels=40), + VacuumField("b_poloidal", "Poloidal Field |B_p|", "T", lower=0.0), + # Signed and unbounded: it diverges wherever B_Z crosses zero, so the + # levels come from a percentile and the stable band is marked with its + # own contours rather than by the colour scale. + VacuumField("decay_index", "Field Decay Index n", "", upper=95.0), + VacuumField("breakdown", "Breakdown Figure E_t B_t / B_p", "V/m", lower=0.0), + ) +} + +#: What ``field=`` may name on the vacuum map. +VACUUM_FIELD_NAMES = tuple(VACUUM_FIELDS) + +#: Where the decay index is passively stable for a rigid current ring. Drawn +#: as two grey contours beneath the filled map, so the band a startup has to +#: sit inside is visible without reading the colourbar. +DECAY_INDEX_STABLE_BAND = (0.0, 1.5) + + +def _vacuum_map_time(ods: Any, time: float | None, time_index: int | None) -> float: + """The instant a vacuum map is drawn at: an index, a time, or the default.""" + if time_index is not None: + base = _array(ods, "pf_active.time") + if base is None or len(base) == 0: + raise ValueError("pf_active.time is required to use time_index=") + index = int(time_index) + if not -len(base) <= index < len(base): + raise ValueError( + f"time_index={time_index} is outside the {len(base)} stored PF " + "samples" + ) + return float(np.asarray(base, dtype=float)[index]) + if time is not None: + return float(time) + return _vacuum_psi_time(ods) + + +def _limiter_interior(ods: Any, r_axis: np.ndarray, z_axis: np.ndarray) -> np.ndarray | None: + """Which ``(R, Z)`` grid points lie inside the limiter, or ``None``. + + A vacuum map's extremes are the coils and the vessel conductors it is + computed from: on VEST several hundred filaments sit inside the limiter's + bounding box, and a grid point landing on one reads that filament's own + singular field. Not one of them lies inside the limiter outline itself, + so the plasma-facing region is exactly the part of the map that means + anything -- and once it is the only part drawn, the contour levels + describe the field a discharge would see instead of a handful of + conductors. + """ + from matplotlib.path import Path as _Path + + layers = _wall_layers(ods) + if not layers: + return None + mesh_r, mesh_z = np.meshgrid(r_axis, z_axis, indexing="ij") + points = np.column_stack([mesh_r.ravel(), mesh_z.ravel()]) + inside = np.zeros(points.shape[0], dtype=bool) + for layer in layers: + outline = np.column_stack([np.asarray(layer.r, dtype=float), + np.asarray(layer.z, dtype=float)]) + if outline.shape[0] >= 3: + inside |= _Path(outline).contains_points(points) + return inside.reshape(mesh_r.shape) if inside.any() else None + + +def _vacuum_levels(values: np.ndarray, field: VacuumField) -> np.ndarray | int: + """Contour levels spanning what the map actually holds. + + ``lower=0`` pins a magnitude's scale to zero, so a null reads as a null + rather than as the bottom of whatever range this instant happens to have. + """ + finite = values[np.isfinite(values)] + if finite.size == 0: + return field.levels + high = float(np.percentile(finite, field.upper)) + low = float(field.lower) if field.lower is not None else float( + np.percentile(finite, 100.0 - field.upper) + ) + if not high > low: + return field.levels + return np.linspace(low, high, field.levels) + + +def _build_vacuum_field( + ods: Any, + *, + field: str = "psi", + time: float | None = None, + time_index: int | None = None, + resolution: int = 65, + units: str | None = None, + **_: Any, +) -> Field2D: + """One quantity of the vacuum field on the poloidal plane, at one instant. + + Flux, poloidal field strength, the decay index and the breakdown figure of + merit are four readings of a single evaluation + (:func:`vaft.omas.process_wrapper.compute_vacuum_field_map`), so they + always describe the same field. The grid's response to the coils and the + vessel is cached, which is what makes ``time_index=`` usable as a slider. + """ + from vaft.formula.equilibrium import ( + decay_index_from_bz, + poloidal_field_magnitude, + toroidal_electric_field, + ) + from vaft.omas.process_wrapper import compute_vacuum_field_map + + if field not in VACUUM_FIELDS: + raise ValueError( + f"unknown vacuum field {field!r}; expected one of " + f"{', '.join(VACUUM_FIELD_NAMES)}" + ) + spec = VACUUM_FIELDS[field] + # Resolve the instant against the caller's whole ODS: the default is the + # breakdown onset, which is read from the magnetics the copy below drops. + instant = _vacuum_map_time(ods, time, time_index) + # The evaluator solves the vessel currents when they are absent; give it a + # private copy so the caller's ODS is left as it was found. + ods = _isolated_copy(ods, (*_NULL_FIELD_ROOTS, "tf")) + result = compute_vacuum_field_map(ods, time=instant, resolution=int(resolution)) + + r_axis, z_axis = result["r"], result["z"] + mesh_r = r_axis[:, None] + if field == "psi": + values = result["psi"] + elif field == "b_poloidal": + values = poloidal_field_magnitude(result["b_r"], result["b_z"]) + elif field == "decay_index": + # Along R, which is axis 0 of an (R, Z) map. + values = decay_index_from_bz(r_axis, result["b_z"], axis=0) + else: + e_toroidal = toroidal_electric_field(mesh_r, result["dpsi_dt"]) + b_toroidal = _vacuum_b_toroidal(ods, result["time"], mesh_r) + b_poloidal = poloidal_field_magnitude(result["b_r"], result["b_z"]) + with np.errstate(divide="ignore", invalid="ignore"): + values = np.abs(e_toroidal) * np.abs(b_toroidal) / b_poloidal + values = np.where(np.isfinite(values), values, np.nan) + + interior = _limiter_interior(ods, r_axis, z_axis) + if interior is not None: + # Masked here rather than through Field2D.region, so the stable-band + # contours are confined too: outside the vessel the decay index is the + # coils' own field and its zero crossings are not a stability boundary. + values = np.where(interior, values, np.nan) + + display = resolve_display( + spec.canonical_unit, unit=units, subject=spec.subject, + quantity=field if not spec.canonical_unit else None, data=values, + ) + scaled = values * display.scale + bracket = f" [{display.unit}]" if display.unit else "" + return Field2D( + # Field2D is laid out (len(z), len(r)); the evaluator answers (R, Z). + r=r_axis, + z=z_axis, + values=scaled.T, + value_label=f"{spec.label}{bracket}", + display=display, + filled=spec.filled, + contour_levels=_vacuum_levels(scaled.T, spec), + secondary_levels=DECAY_INDEX_STABLE_BAND if field == "decay_index" else None, + extend=spec.extend, + overlays=tuple(_wall_layers(ods)), + title=f"{spec.label} at t = {result['time'] * 1e3:.1f} ms (vacuum)", + ) + + +def _vacuum_b_toroidal(ods: Any, time: float, mesh_r: np.ndarray) -> np.ndarray: + """``B_phi = R_0 B_0 / R`` at ``time``, from the TF vacuum field product.""" + product = _array(ods, "tf.b_field_tor_vacuum_r.data") + base = _array(ods, "tf.b_field_tor_vacuum_r.time") + if product is None or base is None or len(product) == 0: + raise ValueError( + "tf.b_field_tor_vacuum_r is required for field='breakdown'; without " + "the toroidal field there is no breakdown figure of merit" + ) + index = int(np.argmin(np.abs(np.asarray(base, dtype=float) - float(time)))) + return float(np.asarray(product, dtype=float)[index]) / mesh_r + + def _build_core_profile_field( ods: Any, *, quantity: str, time_slice: int = 0, **_: Any ) -> Field2D: @@ -3223,6 +3445,11 @@ def _build_coils_non_axisymmetric_topview(ods: Any, **_: Any) -> GeometryLayers: builder=_build_vacuum_psi, description="Vacuum flux map from the PF currents via vaft.omas.compute_null_ods.", ) +RECIPES["vacuum_field"] = CallableRecipe( + builder=_build_vacuum_field, + description="Flux, |B_p|, the decay index or the breakdown figure of merit " + "from one cached vacuum-field evaluation.", +) RECIPES["electron_temperature_field"] = CallableRecipe( builder=lambda ods, **options: _build_core_profile_field( ods, **{**options, "quantity": "temperature"} @@ -5566,6 +5793,10 @@ def _style_psi_field( def field_options_for(name: str) -> tuple[str, ...] | None: """The ``field=`` vocabulary of plot ``name``, or ``None`` if it takes none.""" + if name == "vacuum_field": + # Computed, not read from a stored 2-D array, so it is not a + # FieldRecipe -- but it selects among quantities the same way. + return VACUUM_FIELD_NAMES recipe = RECIPES.get(name) return tuple(recipe.fields) if isinstance(recipe, FieldRecipe) and recipe.fields else None diff --git a/vaft/plot/controls.py b/vaft/plot/controls.py index 4440ef18..54fa9904 100644 --- a/vaft/plot/controls.py +++ b/vaft/plot/controls.py @@ -42,6 +42,10 @@ def _plain(label: str) -> str: #: that is otherwise absent. NONE = "none" +#: Flux maps with no plasma in them: there is no separatrix to normalise +#: against and no magnetic axis to mark, so the psi styles do not apply. +_NO_FLUX_STYLE = frozenset({"equilibrium_field_psi_vacuum", "vacuum_field"}) + @dataclass(frozen=True) class ControlSpec: @@ -142,6 +146,18 @@ def controls_for( def _slice_controls(record: Any) -> list[ControlSpec]: + times: Mapping[str, Any] = getattr(record, "times", None) or {} + if times.get("option") == "time_index" and int(times.get("count", 0)) > 1: + # A dense time base is a slider, not a list: thousands of samples + # cannot be offered as radio buttons, and the reader wants to sweep + # them anyway. The label carries the span, since the positions + # themselves are indices. + count = int(times["count"]) + return [ControlSpec( + "time_index", "range", + f"Time sample ({float(times['start']) * 1e3:.0f}-{float(times['stop']) * 1e3:.0f} ms)", + int(times.get("selected") or 0), (0, count - 1, 1), group="slice", + )] slices: Mapping[str, Any] = getattr(record, "slices", None) or {} usable = tuple(int(i) for i in slices.get("usable", ())) if len(usable) < 2: @@ -236,7 +252,10 @@ def _model_controls(record: Any) -> list[ControlSpec]: default = tuple(name for name in overlay_defaults_for(record.name) if name in overlays) controls.append(ControlSpec("overlay", "multi", "Overlays", default, overlays)) display: Mapping[str, Any] = record.display or {} - if "convention" in display and record.model in ("Field2D", "Panels") and record.name != "equilibrium_field_psi_vacuum": + # The psi styles normalise against the separatrix and mark the axis; a + # vacuum map has neither, so it offers no flux-map style. + if ("convention" in display and record.model in ("Field2D", "Panels") + and record.name not in _NO_FLUX_STYLE): controls.append(ControlSpec( "style", "choice", "Flux map style", PSI_STYLES[0], tuple(PSI_STYLES), applies_to=flux_only, diff --git a/vaft/plot/display.py b/vaft/plot/display.py index 855847c5..ce13bb71 100644 --- a/vaft/plot/display.py +++ b/vaft/plot/display.py @@ -139,7 +139,11 @@ def __post_init__(self) -> None: "current_turns", "A-turns", {"A-turns": 1.0, "kA-turns": 1e-3}, "kA-turns" ), QuantityDisplay("voltage", "V", {"V": 1.0, "mV": 1e3}, "V"), - QuantityDisplay("magnetic_field", "T", {"T": 1.0, "mT": 1e3}, "mT"), + # Gauss is the working unit of tokamak start-up: a breakdown null is a few + # gauss and the vertical field a few hundred, where millitesla reads as a + # decimal fraction. + QuantityDisplay("magnetic_field", "T", {"T": 1.0, "mT": 1e3, "G": 1e4}, "mT"), + QuantityDisplay("electric_field", "V/m", {"V/m": 1.0, "mV/m": 1e3}, "V/m"), QuantityDisplay( "magnetic_flux", "Wb", @@ -207,6 +211,7 @@ def __post_init__(self) -> None: SUBJECT_UNIT_DEFAULTS: dict[tuple[str, str], str] = { ("tf_coil", "magnetic_field"): "T", ("barometry", "pressure"): "Torr", + ("vacuum", "magnetic_field"): "G", } #: (subject, quantity) -> notation override. @@ -222,6 +227,9 @@ def __post_init__(self) -> None: #: is why a beta family plot cannot put all three on one shared axis. DIMENSIONLESS_DISPLAY: dict[tuple[str, str], tuple[str, float, str]] = { ("equilibrium", "beta_t"): ("%", 100.0, "percent"), + # A ratio of a field gradient to the field; the stable window 0 < n < 1.5 + # is quoted in these units and no other, so there is nothing to convert. + ("vacuum", "decay_index"): ("", 1.0, "auto"), } diff --git a/vaft/plot/models.py b/vaft/plot/models.py index 09498b79..b5c5b6da 100644 --- a/vaft/plot/models.py +++ b/vaft/plot/models.py @@ -305,6 +305,12 @@ class Field2D(ViewModel): #: The display policy's resolution of the value unit, when the builder #: applied one; ``value_label`` already carries the unit it names. display: "DisplaySpec | None" = None + #: What happens to values outside ``contour_levels``: ``"neither"`` leaves + #: them blank, ``"min"``/``"max"``/``"both"`` saturate them at the end + #: colours. Levels chosen from a percentile need this -- otherwise the + #: points the percentile deliberately excluded come out as holes in the + #: map, indistinguishable from missing data. + extend: str = "neither" #: Where the main contours are drawn, as a boolean ``(len(z), len(r))`` #: grid; ``None`` draws them everywhere. A flux map confines its plasma #: levels to the plasma, since the same psi values recur beside the coils. @@ -340,6 +346,11 @@ def __post_init__(self) -> None: "contour_levels", as_model_array(self.contour_levels, where="Field2D.contour_levels"), ) + if self.extend not in ("neither", "min", "max", "both"): + raise ValueError( + 'Field2D.extend must be one of "neither", "min", "max", "both"; ' + f"got {self.extend!r}" + ) object.__setattr__(self, "overlays", tuple(self.overlays)) diff --git a/vaft/plot/renderers/fields.py b/vaft/plot/renderers/fields.py index effc91f8..a7c08e9c 100644 --- a/vaft/plot/renderers/fields.py +++ b/vaft/plot/renderers/fields.py @@ -21,6 +21,7 @@ "equilibrium_field_2d", "equilibrium_field_psi", "equilibrium_field_psi_vacuum", + "vacuum_field", "render_field_2d", ] @@ -54,6 +55,8 @@ def render_field_2d( contour_kwargs = {"cmap": cmap, **style} if levels is not None: contour_kwargs["levels"] = levels + if model.extend != "neither": + contour_kwargs["extend"] = model.extend if model.secondary_levels: axes.contour( model.r, model.z, model.values, levels=list(model.secondary_levels), @@ -160,6 +163,27 @@ def equilibrium_field_psi_vacuum( return render_field_2d(model, ax=ax, show=show, **style) +@_field_renderer( + domain="pf_active", quantity="", + subject="vacuum", + description="The vacuum field of the coils and vessel at one instant: the " + "flux, the poloidal field strength, the decay index, or the " + "breakdown figure of merit, chosen with field=.", + ids=("pf_active", "pf_passive", "wall", "tf", "equilibrium", "magnetics"), + required_paths=("pf_active.time", "pf_active.coil.{i}.current.data"), + optional_paths=( + "pf_passive.loop.{i}.element.{j}.geometry.outline.r", + "tf.b_field_tor_vacuum_r.data", + "wall.description_2d.{i}.limiter.unit.{j}.outline.r", + ), +) +def vacuum_field( + model: Field2D, *, ax: Axes | None = None, show: bool = False, **style: Any +) -> tuple[Figure, Axes]: + """One quantity of the coils' and vessel's vacuum field, at one instant.""" + return render_field_2d(model, ax=ax, show=show, **style) + + @_field_renderer( domain="machine", quantity="wall_reduction", subject="passive_structure", diff --git a/vaft/plot/taxonomy.py b/vaft/plot/taxonomy.py index 3945c474..ff80951f 100644 --- a/vaft/plot/taxonomy.py +++ b/vaft/plot/taxonomy.py @@ -102,6 +102,9 @@ class QuantityFamily: Subject("machine", "machine"), # Reconstructions, models, and codes Subject("equilibrium", "reconstruction"), + # The field the coils and the vessel make with no plasma in it: a model of + # the machine, not a reconstruction of a discharge. + Subject("vacuum", "model", ("vacuum_field", "null_field")), Subject("core_profiles", "reconstruction"), Subject("mhd_linear", "model"), Subject("nbi", "machine", ("neutral_beam", "nubeam")), From 8bd059b745496534cec0d2714f1b63cd5e433b74 Mon Sep 17 00:00:00 2001 From: Yun Date: Mon, 7 Sep 2026 18:14:09 +0900 Subject: [PATCH 3/5] Say what the measurement was, not a garbled version of it The level percentile's justification had the ratio backwards: away from a null the map's maximum sits about 30% above its 99th percentile, not within a third of it. Also note that only Matplotlib needs telling about the extend policy -- Plotly clamps to its level range already, so the two backends agree on these maps either way. Co-Authored-By: Claude Opus 5 --- vaft/plot/backend/recipes.py | 12 +++++++----- vaft/plot/models.py | 3 ++- 2 files changed, 9 insertions(+), 6 deletions(-) diff --git a/vaft/plot/backend/recipes.py b/vaft/plot/backend/recipes.py index f32dec8d..c10cde9b 100644 --- a/vaft/plot/backend/recipes.py +++ b/vaft/plot/backend/recipes.py @@ -3054,10 +3054,12 @@ class VacuumField: ``canonical_unit`` is empty for the decay index, which reaches the display policy through :data:`vaft.plot.display.DIMENSIONLESS_DISPLAY` instead. - ``upper`` is the percentile the default contour levels stop at: the field - diverges on the coil filaments the grid sits among, so raw min-to-max - levels would spend the whole colour range on a handful of points and leave - the vessel blank. + ``upper`` is the percentile the default contour levels stop at. Confined + to the limiter these fields are well behaved -- away from a null the map's + maximum sits within about 30% of its 99th percentile -- so stopping there + gives away little, and what falls outside saturates rather than vanishing. + The decay index is the exception: it diverges wherever B_Z crosses zero, + which happens inside the vessel, so it stops sooner. """ name: str @@ -3065,7 +3067,7 @@ class VacuumField: canonical_unit: str subject: str = "vacuum" filled: bool = True - upper: float = 95.0 + upper: float = 99.0 lower: float | None = None levels: int = 24 diff --git a/vaft/plot/models.py b/vaft/plot/models.py index b5c5b6da..c2aed1df 100644 --- a/vaft/plot/models.py +++ b/vaft/plot/models.py @@ -309,7 +309,8 @@ class Field2D(ViewModel): #: them blank, ``"min"``/``"max"``/``"both"`` saturate them at the end #: colours. Levels chosen from a percentile need this -- otherwise the #: points the percentile deliberately excluded come out as holes in the - #: map, indistinguishable from missing data. + #: map, indistinguishable from missing data. Matplotlib needs telling; + #: Plotly clamps to its level range already, so the two agree either way. extend: str = "neither" #: Where the main contours are drawn, as a boolean ``(len(z), len(r))`` #: grid; ``None`` draws them everywhere. A flux map confines its plasma From fc8de35442b53ecddf750cfb979f6970b8681f9f Mon Sep 17 00:00:00 2001 From: Yun Date: Mon, 7 Sep 2026 18:33:53 +0900 Subject: [PATCH 4/5] Declare the IDS the default instant is read from With no time= the map is drawn at the breakdown onset, and that timing reads H-alpha from spectrometer_uv alongside the plasma current. I left it out of the plot's ids, so an adapter that loads only the declared IDSs resolved a different onset than one handed the whole entry -- one PF sample apart, which during the ramp is a 1.3% difference in the drawn flux, not just a different number in the title. Only the IMAS side showed it, because the OMAS adapter passes the whole ODS regardless. equilibrium_field_psi_vacuum has declared both all along, for this exact reason. Found by test_imas_omas_plot_equivalence, and now pinned directly so it cannot come back quietly. The flux equivalence test also has to respect the response budget added after it was written: it now compares against every third point of compute_null_ods's own grid, and bounds a percentile plus the fraction of points that disagree at all rather than a 99th percentile that was really measuring the subsampling. Measured: median 3e-8, 95th percentile 7e-7, 1.4% of points past 1e-3, all of them beside a source filament. Co-Authored-By: Claude Opus 5 --- test/test_vacuum_field_map.py | 25 +++++++++++++++++-------- test/test_vacuum_field_plot.py | 15 +++++++++++++++ vaft/plot/renderers/fields.py | 7 ++++++- 3 files changed, 38 insertions(+), 9 deletions(-) diff --git a/test/test_vacuum_field_map.py b/test/test_vacuum_field_map.py index 140fb497..ef96b05a 100644 --- a/test/test_vacuum_field_map.py +++ b/test/test_vacuum_field_map.py @@ -115,23 +115,32 @@ def test_the_cache_does_not_grow_without_bound(solved): def test_the_flux_matches_compute_null_ods(solved): """Re-pointing the vacuum psi map onto this evaluator must not move it. - The tolerance is a percentile, not a maximum, and deliberately so: the two - paths differ only where a grid point sits on a source filament, where both - answers are meaningless. A max-based bound could not pass and would be - demanding the wrong thing. + The tolerance is a percentile and a count, not a maximum, and deliberately + so: the two paths differ only where a grid point sits on a source filament, + where both answers are meaningless. A max-based bound could not pass and + would be demanding the wrong thing. Measured on the packaged shot: median + 3e-8, 95th percentile 7e-7, and 1.4% of points past 1e-3. """ psi_reference, mesh_r, mesh_z = vaft.omas.compute_null_ods(solved, 0.29) psi_reference = np.asarray(psi_reference, dtype=float) r_axis, z_axis = np.unique(np.asarray(mesh_r)), np.unique(np.asarray(mesh_z)) - result = pw.compute_vacuum_field_map(solved, time=0.29, grid=(r_axis, z_axis)) + # Every third point of the reference grid: the comparison is against its own + # values at its own coordinates, but the whole 129x129 is past the response + # budget -- which is the guard doing its job, not something to work around. + stride = 3 + result = pw.compute_vacuum_field_map( + solved, time=0.29, grid=(r_axis[::stride], z_axis[::stride]) + ) # compute_null_ods lays its grid out (Z, R); this evaluator lays it out (R, Z). - reference = psi_reference.T + reference = psi_reference[::stride, ::stride].T assert reference.shape == result["psi"].shape relative = np.abs(reference - result["psi"]) / np.maximum(np.abs(reference), 1e-12) - assert np.percentile(relative, 95) < 1e-4 - assert np.percentile(relative, 99) < 1e-3 + assert np.percentile(relative, 95) < 1e-5 + # ... and the disagreement is confined to a handful of points, rather than + # being a small bias spread over the map. + assert np.mean(relative > 1e-3) < 0.03 # --------------------------------------------------------------------------- diff --git a/test/test_vacuum_field_plot.py b/test/test_vacuum_field_plot.py index 14bbec96..c5d282ce 100644 --- a/test/test_vacuum_field_plot.py +++ b/test/test_vacuum_field_plot.py @@ -178,3 +178,18 @@ def test_the_caller_s_ods_is_left_as_it_was_found(ods): assert "time" not in ods["pf_passive"] recipes._build_vacuum_field(ods, field="psi", resolution=COARSE) assert "time" not in ods["pf_passive"] + + +def test_the_declared_ids_cover_what_the_default_instant_reads(): + """With no ``time=`` the map is drawn at the breakdown onset, and that + timing reads H-alpha alongside the plasma current. An adapter that loads + only the declared IDSs then resolves a different instant than one handed + the whole entry -- which is a difference in the drawn values, not just in + the title. ``equilibrium_field_psi_vacuum`` declares both for the same + reason. + """ + from vaft.plot.registry import get_spec + + declared = set(get_spec("vacuum_field").ids) + assert {"magnetics", "spectrometer_uv"} <= declared + assert declared >= set(get_spec("equilibrium_field_psi_vacuum").ids) - {"tf"} diff --git a/vaft/plot/renderers/fields.py b/vaft/plot/renderers/fields.py index a7c08e9c..961f50eb 100644 --- a/vaft/plot/renderers/fields.py +++ b/vaft/plot/renderers/fields.py @@ -169,7 +169,12 @@ def equilibrium_field_psi_vacuum( description="The vacuum field of the coils and vessel at one instant: the " "flux, the poloidal field strength, the decay index, or the " "breakdown figure of merit, chosen with field=.", - ids=("pf_active", "pf_passive", "wall", "tf", "equilibrium", "magnetics"), + # spectrometer_uv earns its place: with no time= the map is drawn at the + # breakdown onset, and that timing reads H-alpha alongside the plasma + # current. An adapter that loads only the declared IDSs would otherwise + # resolve a different instant than one that hands over the whole entry. + ids=("pf_active", "pf_passive", "wall", "tf", "equilibrium", "magnetics", + "spectrometer_uv"), required_paths=("pf_active.time", "pf_active.coil.{i}.current.data"), optional_paths=( "pf_passive.loop.{i}.element.{j}.geometry.outline.r", From 265f1892dcf92c5b073a601720cb6e79bbebe69c Mon Sep 17 00:00:00 2001 From: Yun Date: Fri, 11 Sep 2026 15:12:05 +0900 Subject: [PATCH 5/5] Count the three kernels the vacuum map added to the catalog poloidal_field_magnitude, decay_index_from_bz and toroidal_electric_field take the equilibrium category from 107 to 110. Checked against develop's own surface rather than inferred from the failure: exactly those three are new, and nothing develop added was dropped in the rebase. Co-Authored-By: Claude Opus 5 --- test/test_formula_catalog.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/test_formula_catalog.py b/test/test_formula_catalog.py index b389744b..0c0816d6 100644 --- a/test/test_formula_catalog.py +++ b/test/test_formula_catalog.py @@ -57,7 +57,7 @@ def test_the_catalog_counts_the_known_public_surface(): assert counts == { "constants": 0, "utils": 10, - "equilibrium": 107, + "equilibrium": 110, "stability": 19, "green": 16, "atomic": 3,