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, 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..ef96b05a --- /dev/null +++ b/test/test_vacuum_field_map.py @@ -0,0 +1,188 @@ +"""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 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)) + + # 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[::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-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 + + +# --------------------------------------------------------------------------- +# 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/test/test_vacuum_field_plot.py b/test/test_vacuum_field_plot.py new file mode 100644 index 00000000..c5d282ce --- /dev/null +++ b/test/test_vacuum_field_plot.py @@ -0,0 +1,195 @@ +"""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"] + + +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/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/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 110da2ef..bb9672d3 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,265 @@ 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 + +#: 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. + + 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 _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. + + 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) + + +#: 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. + + 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. + _solve_or_recall_eddy_currents(ods) + + 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)))) + + _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] + 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. 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..c10cde9b 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,222 @@ 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. 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 + label: str + canonical_unit: str + subject: str = "vacuum" + filled: bool = True + upper: float = 99.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 +3447,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 +5795,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..c2aed1df 100644 --- a/vaft/plot/models.py +++ b/vaft/plot/models.py @@ -305,6 +305,13 @@ 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. 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 #: levels to the plasma, since the same psi values recur beside the coils. @@ -340,6 +347,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..961f50eb 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,32 @@ 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=.", + # 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", + "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")),