Refactor gk-rz into the naming conventions which are used by GK and G…

.agents/skills/postgkyl-architecture/SKILL.md

@@ -29,11 +29,12 @@ remain GData without importing gdata. Reuse `gdatastate/guards.py` and
 Readers return `(grid, values)` and fill plain metadata; they never import state.
 
 Diagnostics are free functions under model families `gk`, `vm`, `pkpm`, or `mom`,
-not GData methods. Model-specific loading belongs beside its physics. Resolve
+not GData methods. Model-specific loading and auxiliary-file discovery belong
+beside their compositions. Coordinate operations receive explicit geometry or projections;
+operations has no equation-family subpackages. Resolve
 output stems/frames through `diagnostics/discovery.py`; keep quantity vocabulary
 in the equation module's `VARIABLES` table. Use public functions from lower
 layers and return a state via `_result`, or a Figure for program diagnostics.
-for the current major version.
 
 Run the import, foreign-floor, facade, and canonical-callable contracts in
 `tests/test_postgkyl.py` after structural changes.

docs/source/physical-rz.rst

@@ -9,6 +9,13 @@ The Python and CLI outputs below are both generated and compared when this websi
 Computational coordinates are useful for analysis, but a poloidal view places
 the field in its physical geometry. Keep the companion
 ``rt_gk_tcv_nt_iwl_3x2v_p1-geo_int_mapc2p.gkyl`` beside the electron density
-file so the operation can resolve it by the simulation prefix.
+file so ``pg.gk.rz`` can resolve it by the simulation prefix. This diagnostic
+loads Gkeyll geometry and composes equation-independent coordinate operations.
+
+For explicit geometry, construct ``pg.Geometry`` from coordinate arrays, build
+``pg.resolve_rz_projection(data, geometry)``, and apply it with
+``data.map_to_rz(projection=projection)``. The projection can be reused for
+fields on the same computational grid. Three-dimensional reconstruction
+assumes periodic field-aligned coordinates and twist-and-shift boundaries.
 
 .. include:: _pairs/05_gk_rz.inc

examples/README.md

@@ -20,7 +20,7 @@ of truth for "does the tutorial still work," not just this README.
 | [`02_arithmetic_and_numpy.py`](scripts/02_arithmetic_and_numpy.py) | Weak DG algebra on raw modal data (`*`, `/`, `+`, `.integrate()`) vs. plain NumPy math after `.interpolate()`, and the guardrail between them |
 | [`03_diagnostics_five_moment.py`](scripts/03_diagnostics_five_moment.py) | The `diagnostics` layer: equation-specific physics (`postgkyl.diagnostics.mom.five_moment`) on top of a `GData`, on a generated shock-tube initial state |
 | [`04_gyrokinetics.py`](scripts/04_gyrokinetics.py) | The gyrokinetic diagnostics: `pg.gk.load_quantity` (named moments/geometry, resolved by naming convention) and `pg.gk.load_distf` (full distribution function), on the `rt_gk_tcv_iwl*` fixtures |
-| [`05_gk_rz.py`](scripts/05_gk_rz.py) | The gyrokinetic R-Z operation: one-line fluent and functional calls plus projection reuse over multiple toroidal angles |
+| [`05_gk_rz.py`](scripts/05_gk_rz.py) | Gkeyll geometry discovery with `pg.gk.rz`, plus explicit projection reuse with `map_to_rz` |
 | [`06_growth.py`](scripts/06_growth.py) | Recover a known energy growth rate and inspect log-space residuals |
 | [`07_collect_animate.py`](scripts/07_collect_animate.py) | Load many frames, collect a space–time diagram, and animate 1D/2D travelling waves |
 | [`08_plotly.py`](scripts/08_plotly.py) | Interactive 2D height surfaces and 3D volume isosurfaces |
@@ -49,7 +49,7 @@ See [`cli_tutorial.md`](cli_tutorial.md) -- inspecting a file, the
 `interpolate`/`select`/`plot` chain, discontinuity-preserving plots with
 `local_poly`, DynVector `info`/`fit`, the gyrokinetic loaders
 (`gk_load_quantity`, `gk_load_distf`), the `gk_rz`
-transformation, `save`, and the generated command inventory.
+diagnostic, `save`, and the generated command inventory.
 
 ## Running the tests
 

examples/cli_tutorial.md

@@ -125,10 +125,10 @@ pgkyl gk_load_distf --name tests/test_data/rt_gk_tcv_iwl_1x2v_p1 \
 
 ## 7. Map a gyrokinetic field to R-Z
 
-`gk_rz` is a data transformation: it interpolates one raw DG component and
+`gk_rz` is a diagnostic composition: it interpolates one raw DG component and
 maps it onto the physical poloidal plane. Geometry is inferred from the
 field's filename, preferring nodal geometry and falling back to modal
-`mapc2p` geometry. The CLI and Python calls below use the same operation and
+`mapc2p` geometry. The CLI and Python calls below use the same diagnostic and
 defaults:
 
 ```bash
@@ -139,9 +139,10 @@ pgkyl tests/test_data/rt_gk_tcv_nt_iwl_3x2v_p1-elc_M0_5.gkyl \
 ```python
 import postgkyl as pg
 
-mapped = pg.load(
+data = pg.load(
     "tests/test_data/rt_gk_tcv_nt_iwl_3x2v_p1-elc_M0_5.gkyl"
-).gk_rz(nz_interp=2)
+)
+mapped = pg.gk.rz(data, nz_interp=2)
 ```
 
 ## 8. Saving to another format

examples/scripts/05_gk_rz.py

@@ -12,7 +12,6 @@
 import numpy as np
 
 import postgkyl as pg
-from postgkyl.operations import gyrokinetics as gk_ops
 
 from _example_paths import TEST_DATA, prepare_output_dir
 
@@ -21,7 +20,7 @@
 data = pg.load(TEST_DATA / "rt_gk_tcv_nt_iwl_3x2v_p1-elc_M0_5.gkyl")
 
 # The common path: geometry is inferred from FIELD's simulation prefix.
-mapped = data.gk_rz(z_axis=0.0, phi_tor=0.0, nz_interp=2)
+mapped = pg.gk.rz(data, z_axis=0.0, phi_tor=0.0, nz_interp=2)
 fig = mapped.plot(title="Electron density in the poloidal plane",
                   xlabel="R [m]",
                   ylabel="Z [m]",
@@ -30,19 +29,12 @@
                   no_show=True)
 fig.savefig(OUTPUT_DIR / "05_gk_rz.png")
 
-# The functional spelling is the identical operation.
-functional = pg.gk_rz(data, z_axis=0.0, phi_tor=0.0, nz_interp=2)
-np.testing.assert_allclose(functional.values, mapped.values)
-
 # Reuse geometry and projection when several fields/frames or toroidal angles
 # share one computational grid.
-geometry = gk_ops.resolve_geometry(data.file_name)
-projection = gk_ops.resolve_rz_projection(data,
-                                          geometry,
-                                          z_axis=0.0,
-                                          nz_interp=2)
-at_zero = gk_ops.map_to_rz(data, projection, phi_tor=0.0)
-at_quarter_turn = gk_ops.map_to_rz(data, projection, phi_tor=np.pi / 2)
+geometry = pg.gk.resolve_geometry(data.file_name)
+projection = pg.resolve_rz_projection(data, geometry, z_axis=0.0, nz_interp=2)
+at_zero = data.map_to_rz(projection=projection, phi_tor=0.0)
+at_quarter_turn = pg.map_to_rz(data, projection=projection, phi_tor=np.pi / 2)
 np.testing.assert_allclose(at_zero.values, mapped.values)
 assert not np.allclose(
     at_zero.values, at_quarter_turn.values, rtol=1e-12, atol=0.0)

src/postgkyl/__init__.py

@@ -17,7 +17,7 @@
     info                             <- operations/ (the info verb, one-or-many)
     integrate                        <- operations/ (grid integral, via Gkeyll)
     interpolate, select              <- operations/ (functional verb spellings)
-    gk_rz                           <- operations/gyrokinetics/ (domain operation)
+    map_to_rz                       <- operations/map.py
     represent, apply                 <- operations/ (value_form verbs)
     available_evaluate_operators     <- operations/ (``evaluate``'s RPN token vocabulary)
     save                             <- io/        (file output)
@@ -75,7 +75,10 @@
     select,
     val2coord,
 )
-from postgkyl.operations.gyrokinetics import gk_fluxsurf, gk_rz
+from postgkyl.operations import (Geometry, RzProjection, FluxSurfaceGrid,
+                                 map_to_rz, resolve_rz_projection,
+                                 extract_flux_surface,
+                                 resolve_flux_surface_grid)
 from postgkyl.render import animate, plot, plotly, plotly_animate, pyvista
 from postgkyl.gdatastate import group_blocks
 from postgkyl.cli_spec import hidden
@@ -93,7 +96,9 @@
     "integrate", "interpolate", "local_poly", "select", "average",
     "eval_at_coord_proj", "fft", "magsq", "mask", "grid", "val2coord",
     "extract_input", "fit", "growth", "differentiate", "map", "represent",
-    "apply", "gk_rz", "gk_fluxsurf", "save", "collect", "evaluate", "relchange",
+    "apply", "Geometry", "RzProjection", "FluxSurfaceGrid", "map_to_rz",
+    "resolve_rz_projection", "extract_flux_surface",
+    "resolve_flux_surface_grid", "save", "collect", "evaluate", "relchange",
     "animate", "plotly_animate", "sort", "available_evaluate_operators",
     "plotly", "pyvista", "gk", "__version__", "version_report"
 ]

src/postgkyl/cli/discovery.py

@@ -102,10 +102,7 @@ def consider(obj, path: str, name: str) -> None:
     namespace = relative.rsplit(".", 1)[-1]
     for name, value in _functions(module):
       if not value.__module__.startswith("postgkyl.diagnostics"):
-        # Compatibility re-exports owned by a lower layer keep that owner's
-        # canonical command (for example operations.gyrokinetics.gk_rz ->
-        # ``gk_rz``); the diagnostic alias is classified but does not invent
-        # a second command or move it into the wrong help section.
+        # Lower-layer re-exports retain their canonical command vocabulary.
         _classify(value, f"{module.__name__}.{name}")
         continue
       consider(value, f"{module.__name__}.{name}", f"{namespace}_{name}")

src/postgkyl/diagnostics/__init__.py

@@ -5,9 +5,8 @@
 here take loaded ``GData``/``GDataState`` (one or several) plus physical
 scalars as keyword-only options, and return a ``GDataState`` (via
 ``_result``) or, in later layers, a ``Figure``. Equation-blind core verbs
-stay in flat ``operations`` modules. Domain-specific transformations live in
-operation subpackages (for example ``operations.gyrokinetics``); this layer
-is reserved for code that knows what field components physically mean.
+stay in flat ``operations`` modules. Model-specific auxiliary discovery and
+interpretation live beside their diagnostic compositions.
 
 The four public packages mirror Gkeyll's model families: ``gk``, ``vm``,
 ``pkpm``, and ``mom``. The equation-blind ``discovery`` module
@@ -150,11 +149,7 @@ def _combine(function, dataset_names: tuple[str, ...]) -> None:
     "is_geo_mapc2p",
     "multib_tag",
     "nodes_to_RZ",
-    "map_to_rz",
     "resolve_geometry",
-    "resolve_rz_projection",
-    "extract_flux_surface",
-    "resolve_flux_surface_grid",
 ):
   _function = getattr(gk, _name)
   if command_spec(_function) is None and hidden_spec(_function) is None:

src/postgkyl/diagnostics/gk/__init__.py

@@ -7,9 +7,8 @@
 instruction file's decision record): splitting resolution from physics would
 give gyrokinetics two homes for one piece of equation knowledge. Only the
 equation-blind stem/frame discovery is shared, via
-``postgkyl.diagnostics.discovery``. Geometry-only transformations live below
-this physics layer in ``postgkyl.operations.gyrokinetics``; the R-Z and
-flux-surface names exported here are compatibility aliases.
+``postgkyl.diagnostics.discovery``. Geometry discovery lives here and composes
+explicit coordinate transformations from ``operations.map``.
 """
 
 from __future__ import annotations
@@ -37,12 +36,11 @@
 # Layer 13: program-scale diagnostics ported from src_bak's apps/gk_*.py.
 from .energy_balance import EnergyBalanceTraces, energy_balance_error, energy_balance
 from .particle_balance import ParticleBalanceTraces, particle_balance, particle_balance_error
-from .nodes import GKYL_GEOMETRY_ID, nodes, is_geo_mapc2p, multib_tag, nodes_to_RZ
+from .nodes import nodes, multib_tag, nodes_to_RZ
 
-# Compatibility exports: canonical transformation APIs now live under
-# postgkyl.operations.gyrokinetics. These imports are exact aliases.
-from .rz import Geometry, RzProjection, gk_rz, map_to_rz, resolve_geometry, resolve_rz_projection
-from .fluxsurf import FluxSurfaceGrid, extract_flux_surface, resolve_flux_surface_grid
+from .geometry import GKYL_GEOMETRY_ID, is_geo_mapc2p, resolve_geometry
+from .rz import rz
+from .fluxsurf import fluxsurf
 
 from typing import Annotated
 
@@ -69,6 +67,9 @@
 for _function in (energy_balance, particle_balance, nodes):
   command(_REPORT_SPEC)(_function)
 
+for _function in (rz, fluxsurf):
+  command(CommandSpec(Section.DIAGNOSTICS, Execution.MAP_REPLACE))(_function)
+
 __all__ = [
     "load_distf",
     "resolve_frames",
@@ -100,13 +101,7 @@
     "is_geo_mapc2p",
     "multib_tag",
     "nodes_to_RZ",
-    "Geometry",
-    "RzProjection",
-    "gk_rz",
-    "map_to_rz",
     "resolve_geometry",
-    "resolve_rz_projection",
-    "FluxSurfaceGrid",
-    "extract_flux_surface",
-    "resolve_flux_surface_grid",
+    "rz",
+    "fluxsurf",
 ]

src/postgkyl/diagnostics/gk/fluxsurf.py

@@ -1,23 +1,83 @@
-"""Compatibility aliases for gyrokinetic flux-surface operations.
-
-Canonical imports live in :mod:`postgkyl.operations.gyrokinetics`; this path
-is scheduled for removal in the next major version.
-"""
-
-from postgkyl.operations.gyrokinetics.fluxsurf import (
-    FluxSurfaceGrid,
-    Geometry,
-    extract_flux_surface,
-    flux_surface_grids,
-    grid_for,
-    resolve_flux_surface_grid,
-)
-
-__all__ = [
-    "Geometry",
-    "FluxSurfaceGrid",
-    "extract_flux_surface",
-    "flux_surface_grids",
-    "grid_for",
-    "resolve_flux_surface_grid",
-]
+"""Gkeyll geometry discovery composed with explicit fluxsurf mapping."""
+from __future__ import annotations
+
+from postgkyl.gdatastate.gdatastate import GDataState
+from postgkyl.operations.map import resolve_flux_surface_grid, extract_flux_surface
+from postgkyl.operations.geometry import FluxSurfaceGrid, validate_mapping_grid
+from .geometry import geometry_prefix, per_block_path, resolve_geometry
+
+
+def flux_surface_grids(datasets,
+                       *,
+                       mapc2p: str | None = None,
+                       nodes_file: str | None = None,
+                       x_idx: int = 0,
+                       nphi: int = 128,
+                       nz_interp: int = 8) -> dict[str | None, FluxSurfaceGrid]:
+  """Build one reusable flux-surface grid per block geometry."""
+  grids: dict[str | None, FluxSurfaceGrid] = {}
+  for data in datasets:
+    key = geometry_prefix(data.file_name)
+    if key in grids:
+      validate_mapping_grid(data, grids[key].computational_grid)
+      continue
+    block = data.ctx.get("block")
+    geometry = resolve_geometry(data.file_name,
+                                mapc2p=per_block_path(mapc2p, block),
+                                nodes_file=per_block_path(nodes_file, block))
+    grids[key] = resolve_flux_surface_grid(data,
+                                           geometry,
+                                           x_idx=x_idx,
+                                           nphi=nphi,
+                                           nz_interp=nz_interp)
+  return grids
+
+
+def grid_for(grids: dict[str | None, FluxSurfaceGrid],
+             data: "GDataState") -> FluxSurfaceGrid:
+  """Return the sampling grid belonging to ``data``'s block."""
+  return grids[geometry_prefix(data.file_name)]
+
+
+def fluxsurf(data: "GDataState",
+             *,
+             mapc2p: str | None = None,
+             nodes_file: str | None = None,
+             x_idx: int = 0,
+             nphi: int = 128,
+             nz_interp: int = 8,
+             comp: int = 0,
+             inplace: bool = False,
+             tag: str | None = None,
+             label: str | None = None) -> "GDataState":
+  """Extract one field component on a toroidal flux surface.
+
+  Args:
+    data: Three-dimensional field-aligned modal dataset.
+    mapc2p: Explicit modal geometry path.
+    nodes_file: Explicit nodal geometry path.
+    x_idx: Radial cell index identifying the surface.
+    nphi: Number of toroidal-angle slices.
+    nz_interp: Parallel-direction interpolation factor.
+    comp: Physical field component to extract.
+    inplace: Mutate and return ``data`` instead of creating a dataset.
+    tag: Optional tag for the returned dataset.
+    label: Optional label for the returned dataset.
+  """
+  geometry = resolve_geometry(data.file_name,
+                              mapc2p=mapc2p,
+                              nodes_file=nodes_file)
+  fs_grid = resolve_flux_surface_grid(data,
+                                      geometry,
+                                      x_idx=x_idx,
+                                      nphi=nphi,
+                                      nz_interp=nz_interp)
+  return extract_flux_surface(data,
+                              fs_grid=fs_grid,
+                              comp=comp,
+                              inplace=inplace,
+                              tag=tag,
+                              label=label)
+
+
+__all__ = ["fluxsurf", "flux_surface_grids", "grid_for"]