From db9e200158596c41bd5afce24a1c88f3e935634a Mon Sep 17 00:00:00 2001
From: munechika-koyo <munechika.koyo@gmail.com>
Date: Mon, 3 Aug 2026 16:37:38 +0200
Subject: [PATCH 01/36] =?UTF-8?q?=E2=9C=A8=20Refactor=20load=5Fgrid=20func?=
 =?UTF-8?q?tion=20to=20support=20entry=20parameter=20and=20improve=20grid?=
 =?UTF-8?q?=20source=20resolution;=20add=20path=20helper=20functions=20for?=
 =?UTF-8?q?=20IDS=20references?=
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit

---
 src/cherab/imas/ids/common/__init__.py      |  10 +
 src/cherab/imas/ids/common/ggd/load_grid.py |  98 +++++++--
 src/cherab/imas/ids/common/path.py          | 207 ++++++++++++++++++++
 3 files changed, 303 insertions(+), 12 deletions(-)
 create mode 100644 src/cherab/imas/ids/common/path.py

diff --git a/src/cherab/imas/ids/common/__init__.py b/src/cherab/imas/ids/common/__init__.py
index db058e1..17e2b0b 100644
--- a/src/cherab/imas/ids/common/__init__.py
+++ b/src/cherab/imas/ids/common/__init__.py
@@ -22,6 +22,12 @@
 from . import ggd, grid_radial, species
 from ._ids_numeric import get_ids_numeric_field
 from ._model import solve_coronal_equilibrium
+from .path import (
+    IDSPathReference,
+    load_ids_path_reference,
+    parse_ids_path_fragment,
+    resolve_ids_path_reference,
+)
 from .slice import get_ids_time_slice
 
 __all__ = [
@@ -31,4 +37,8 @@
     "solve_coronal_equilibrium",
     "get_ids_time_slice",
     "get_ids_numeric_field",
+    "IDSPathReference",
+    "parse_ids_path_fragment",
+    "resolve_ids_path_reference",
+    "load_ids_path_reference",
 ]
diff --git a/src/cherab/imas/ids/common/ggd/load_grid.py b/src/cherab/imas/ids/common/ggd/load_grid.py
index ab29726..ad27b87 100644
--- a/src/cherab/imas/ids/common/ggd/load_grid.py
+++ b/src/cherab/imas/ids/common/ggd/load_grid.py
@@ -17,15 +17,21 @@
 # under the Licence.
 """Module for loading GGD grids from IMAS grid_ggd IDS structures."""
 
-from typing import Literal, overload
+from __future__ import annotations
+
+from typing import Literal, cast, overload
 
 from numpy import int32
 from numpy.typing import NDArray
 
+from imas.db_entry import DBEntry
+from imas.ids_struct_array import IDSStructArray
 from imas.ids_structure import IDSStructure
+from imas.ids_toplevel import IDSToplevel
 
 from ....ggd.unstruct_2d_extend_mesh import UnstructGrid2DExtended
 from ....ggd.unstruct_2d_mesh import UnstructGrid2D
+from ..path import load_ids_path_reference, resolve_ids_path_reference
 from .load_unstruct_2d import load_unstruct_grid_2d
 from .load_unstruct_3d import load_unstruct_grid_2d_extended
 
@@ -34,18 +40,30 @@
 
 @overload
 def load_grid(
-    grid_ggd: IDSStructure, with_subsets: Literal[False] = False, num_toroidal: int | None = None
+    grid_ggd: IDSStructure,
+    with_subsets: Literal[False] = False,
+    num_toroidal: int | None = None,
+    *,
+    entry: DBEntry | None = None,
 ) -> UnstructGrid2D | UnstructGrid2DExtended: ...
 
 
 @overload
 def load_grid(
-    grid_ggd: IDSStructure, with_subsets: Literal[True], num_toroidal: int | None = None
+    grid_ggd: IDSStructure,
+    with_subsets: Literal[True],
+    num_toroidal: int | None = None,
+    *,
+    entry: DBEntry | None = None,
 ) -> tuple[UnstructGrid2D, dict[str, NDArray[int32]], dict[str, int]]: ...
 
 
 def load_grid(
-    grid_ggd: IDSStructure, with_subsets: bool = False, num_toroidal: int | None = None
+    grid_ggd: IDSStructure,
+    with_subsets: bool = False,
+    num_toroidal: int | None = None,
+    *,
+    entry: DBEntry | None = None,
 ) -> (
     UnstructGrid2D
     | tuple[UnstructGrid2D, dict[str, NDArray[int32]], dict[str, int]]
@@ -69,6 +87,8 @@ def load_grid(
         Read grid subset data, by default is False.
     num_toroidal
         Number of toroidal points, by default None.
+    entry
+        Open IMAS data entry used to resolve ``grid_ggd.path`` references to external IDSs.
 
     Returns
     -------
@@ -100,14 +120,15 @@ def load_grid(
 
         grid, subsets, subset_id = load_grid(ids.grid_ggd[0], with_subsets=True)
     """
-    spaces = get_standard_spaces(grid_ggd)
+    grid_source = _resolve_grid_source(grid_ggd, entry=entry)
+    spaces = get_standard_spaces(grid_source)
 
     if not len(spaces):
         raise RuntimeError("GGD grid contain no spaces.")
 
     if len(spaces) == 1:  # simple unstructured grids
         if len(spaces[0].objects_per_dimension) == 3:  # 2D case
-            return load_unstruct_grid_2d(grid_ggd, 0, with_subsets=with_subsets)
+            return load_unstruct_grid_2d(grid_source, 0, with_subsets=with_subsets)
         if len(spaces[0].objects_per_dimension) == 4:  # 3D case
             raise NotImplementedError(
                 "Loading unstructured 3D grids will be implemented in the future."
@@ -118,7 +139,7 @@ def load_grid(
     if len(spaces) == 2:  # 2D structured grid or 2D unstructured grid extended in 3D
         if len(spaces[0].objects_per_dimension) == 3 and len(spaces[1].objects_per_dimension) < 3:
             return load_unstruct_grid_2d_extended(
-                grid_ggd, with_subsets=with_subsets, num_toroidal=num_toroidal
+                grid_source, with_subsets=with_subsets, num_toroidal=num_toroidal
             )
         if len(spaces[0].objects_per_dimension) < 3 and len(spaces[1].objects_per_dimension) < 3:
             raise NotImplementedError(
@@ -142,6 +163,62 @@ def load_grid(
     raise RuntimeError("Unsupported grid type.")
 
 
+def _resolve_grid_source(grid_ggd: IDSStructure, entry: DBEntry | None = None) -> IDSStructure:
+    if len(grid_ggd.space):
+        return grid_ggd
+
+    if not len(grid_ggd.path):
+        return grid_ggd
+
+    path = str(grid_ggd.path).strip()
+    if not path:
+        return grid_ggd
+
+    resolved: IDSToplevel | IDSStructure | IDSStructArray
+
+    if "#" in path:
+        if entry is None:
+            raise RuntimeError(
+                "Unable to resolve grid_ggd.path without a DBEntry. "
+                + "Pass an open entry via load_grid(..., entry=entry)."
+            )
+        resolved = load_ids_path_reference(entry, path)
+    else:
+        root = _find_toplevel(grid_ggd)
+        ids_name = root.metadata.name
+        if path.startswith("/"):
+            resolved = resolve_ids_path_reference(root, f"#{ids_name}{path}")
+        else:
+            resolved = resolve_ids_path_reference(root, path)
+
+    if isinstance(resolved, IDSStructArray):
+        if not len(resolved):
+            raise RuntimeError(f"Resolved grid reference '{path}' points to an empty array.")
+        resolved = resolved[0]
+
+    if not hasattr(resolved, "space"):
+        raise RuntimeError(
+            f"Resolved grid reference '{path}' does not point to a grid_ggd structure."
+        )
+
+    return cast(IDSStructure, resolved)
+
+
+def _find_toplevel(node: IDSStructure) -> IDSToplevel:
+    current = node
+    while hasattr(current, "_parent"):
+        parent = current._parent
+        if isinstance(parent, IDSToplevel):
+            return parent
+        if not hasattr(parent, "_parent"):
+            break
+        current = parent
+
+    raise RuntimeError(
+        "Unable to resolve local grid_ggd.path: cannot find the parent IDS toplevel object."
+    )
+
+
 def get_standard_spaces(grid_ggd: IDSStructure) -> list[IDSStructure]:
     """Get a list of standard non-empty spaces from the grid_ggd structure.
 
@@ -157,14 +234,11 @@ def get_standard_spaces(grid_ggd: IDSStructure) -> list[IDSStructure]:
 
     Raises
     ------
-    ValueError
+    RuntimeError
         If no spaces are defined in the grid_ggd structure.
     """
     if not len(grid_ggd.space):
-        error_massage = "Unable to read the grid. Grid space is not defined."
-        if len(grid_ggd.path):
-            error_massage += f" The grid is defined in {grid_ggd.path}."
-        raise ValueError(error_massage)
+        raise RuntimeError("Unable to read the grid. Grid space is not defined.")
 
     # Get list of standard spaces:
     spaces = []
diff --git a/src/cherab/imas/ids/common/path.py b/src/cherab/imas/ids/common/path.py
new file mode 100644
index 0000000..f0da271
--- /dev/null
+++ b/src/cherab/imas/ids/common/path.py
@@ -0,0 +1,207 @@
+"""Helpers for parsing and resolving IMAS IDS path references."""
+
+from __future__ import annotations
+
+from dataclasses import dataclass
+from typing import cast
+
+from imas.db_entry import DBEntry
+from imas.ids_struct_array import IDSStructArray
+from imas.ids_structure import IDSStructure
+from imas.ids_toplevel import IDSToplevel
+
+__all__ = [
+    "IDSPathReference",
+    "parse_ids_path_fragment",
+    "resolve_ids_path_reference",
+    "load_ids_path_reference",
+]
+
+
+@dataclass(frozen=True, slots=True)
+class IDSPathReference:
+    """Representation of an IMAS same-document IDS path reference."""
+
+    ids_name: str
+    """Referenced IDS name."""
+
+    occurrence: int | None = None
+    """Referenced IDS occurrence, if explicitly provided."""
+
+    idspath: str = ""
+    """IDS path fragment that follows the referenced IDS name."""
+
+
+def parse_ids_path_fragment(reference: str) -> IDSPathReference:
+    """Parse a same-document IMAS IDS path fragment.
+
+    Parameters
+    ----------
+    reference
+        A string in the form ``#ids[:occurrence][/idspath]``.
+
+    Returns
+    -------
+    `.IDSPathReference`
+        Parsed IDS path reference.
+
+    Raises
+    ------
+    ValueError
+        If the reference is empty or does not follow the expected fragment syntax.
+    """
+    if not reference:
+        raise ValueError("The IDS path reference cannot be empty.")
+
+    fragment = reference.rsplit("#", 1)[-1].strip()
+    if not fragment:
+        raise ValueError("The IDS path reference cannot be empty.")
+
+    ids_part, separator, idspath = fragment.partition("/")
+
+    if not ids_part:
+        raise ValueError(f"Invalid IDS path reference '{reference}'.")
+
+    ids_name, colon, occurrence_text = ids_part.partition(":")
+    if not ids_name:
+        raise ValueError(f"Invalid IDS path reference '{reference}'.")
+
+    occurrence: int | None
+    if colon:
+        if not occurrence_text:
+            raise ValueError(f"Invalid IDS path reference '{reference}'.")
+        try:
+            occurrence = int(occurrence_text)
+        except ValueError as exc:
+            raise ValueError(f"Invalid IDS occurrence '{occurrence_text}'.") from exc
+    else:
+        occurrence = None
+
+    return IDSPathReference(
+        ids_name=ids_name,
+        occurrence=occurrence,
+        idspath=f"/{idspath}" if separator else "",
+    )
+
+
+def resolve_ids_path_reference(
+    ids: IDSToplevel | IDSStructure | IDSStructArray,
+    reference: IDSPathReference | str,
+) -> IDSToplevel | IDSStructure | IDSStructArray:
+    """Resolve an IDS path reference against an already loaded IDS object.
+
+    Parameters
+    ----------
+    ids
+        The loaded IDS root object to resolve the path against.
+    reference
+        Parsed IDS path reference or a string in the form ``#ids[:occurrence][/idspath]``.
+
+    Returns
+    -------
+    `~imas.ids_toplevel.IDSToplevel` | `~imas.ids_structure.IDSStructure` | `~imas.ids_struct_array.IDSStructArray`
+        The resolved IDS object.
+
+    Raises
+    ------
+    ValueError
+        If the reference does not target the provided IDS root, or the path syntax is invalid.
+    """
+    ref = parse_ids_path_fragment(reference) if isinstance(reference, str) else reference
+
+    ids_name = _get_ids_name(ids)
+    if ids_name is not None and ref.ids_name != ids_name:
+        raise ValueError(
+            f"The reference targets IDS '{ref.ids_name}', but the provided IDS root is '{ids_name}'."
+        )
+
+    if not ref.idspath:
+        return ids
+
+    return _resolve_idspath(ids, ref.idspath)
+
+
+def load_ids_path_reference(
+    entry: DBEntry,
+    reference: IDSPathReference | str,
+) -> IDSToplevel | IDSStructure | IDSStructArray:
+    """Load and resolve an IDS path reference from a data entry.
+
+    Parameters
+    ----------
+    entry
+        Open IMAS data entry.
+    reference
+        Parsed IDS path reference or a string in the form ``#ids[:occurrence][/idspath]``.
+
+    Returns
+    -------
+    `~imas.ids_toplevel.IDSToplevel` | `~imas.ids_structure.IDSStructure` | `~imas.ids_struct_array.IDSStructArray`
+        The resolved IDS object.
+    """
+    ref = parse_ids_path_fragment(reference) if isinstance(reference, str) else reference
+    ids = entry.get(ref.ids_name, occurrence=ref.occurrence or 0, autoconvert=False, lazy=True)
+    if not ref.idspath:
+        return ids
+
+    return _resolve_idspath(ids, ref.idspath)
+
+
+def _get_ids_name(ids: IDSToplevel | IDSStructure | IDSStructArray) -> str | None:
+    metadata = getattr(ids, "metadata", None)
+    return getattr(metadata, "name", None)
+
+
+def _resolve_idspath(
+    ids: IDSToplevel | IDSStructure | IDSStructArray,
+    idspath: str,
+) -> IDSToplevel | IDSStructure | IDSStructArray:
+    current: IDSToplevel | IDSStructure | IDSStructArray = ids
+
+    for segment in idspath.lstrip("/").split("/"):
+        if not segment:
+            raise ValueError(f"Invalid IDS path '{idspath}'.")
+
+        name, index_expression = _split_segment(segment)
+        current = cast(IDSToplevel | IDSStructure | IDSStructArray, getattr(current, name))
+
+        if index_expression is None:
+            continue
+
+        if not isinstance(current, IDSStructArray):
+            raise ValueError(f"Segment '{segment}' refers to a non-array IDS node.")
+
+        current = _resolve_struct_array_index(current, index_expression)
+
+    return current
+
+
+def _split_segment(segment: str) -> tuple[str, str | None]:
+    if "(" not in segment:
+        return segment, None
+
+    name, _, remainder = segment.partition("(")
+    if not name or not remainder.endswith(")"):
+        raise ValueError(f"Invalid IDS path segment '{segment}'.")
+
+    return name, remainder[:-1]
+
+
+def _resolve_struct_array_index(
+    array: IDSStructArray, index_expression: str
+) -> IDSStructure | IDSStructArray:
+    index_expression = index_expression.strip()
+
+    if not index_expression or index_expression == ":":
+        return array
+
+    if index_expression.startswith("{") or ":" in index_expression:
+        raise NotImplementedError(
+            "Resolving IDS path fragments with array slices or index sets is not supported yet."
+        )
+
+    index = int(index_expression)
+    if index == 0:
+        raise ValueError("IMAS IDS path indices are 1-based and cannot be 0.")
+
+    return array[index - 1 if index > 0 else index]

From 81c58b7ca548b6a050f358e084bff70658146db3 Mon Sep 17 00:00:00 2001
From: munechika-koyo <munechika.koyo@gmail.com>
Date: Mon, 3 Aug 2026 16:38:00 +0200
Subject: [PATCH 02/36] =?UTF-8?q?=F0=9F=8E=A8=20Enhance=20grid=20loading?=
 =?UTF-8?q?=20and=20entry=20reference=20handling=20in=20plasma=20modules;?=
 =?UTF-8?q?=20add=20get=5Fentry=5Freference=20utility=20function?=
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit

---
 src/cherab/imas/emitter/radiation.py | 83 ++++++++++++++++++++++++----
 src/cherab/imas/plasma/blend.py      | 14 ++++-
 src/cherab/imas/plasma/core.py       | 10 +++-
 src/cherab/imas/plasma/edge.py       | 26 +++++----
 src/cherab/imas/plasma/utility.py    | 31 +++++++++++
 5 files changed, 136 insertions(+), 28 deletions(-)

diff --git a/src/cherab/imas/emitter/radiation.py b/src/cherab/imas/emitter/radiation.py
index e3791ea..2177818 100644
--- a/src/cherab/imas/emitter/radiation.py
+++ b/src/cherab/imas/emitter/radiation.py
@@ -41,14 +41,14 @@
 
 from ..ggd import UnstructGrid2DExtended
 from ..ggd.base_mesh import InterpolatorCacheMode
-from ..ids.common import get_ids_time_slice
+from ..ids.common import get_ids_time_slice, load_ids_path_reference, resolve_ids_path_reference
 from ..ids.common.ggd import load_grid
 from ..ids.common.grid_radial import GridData, get_psi_norm, load_core_grid
 from ..ids.radiation import load_core_emissivity, load_ggd_emissivity
 from ..math import FourierBezierConstructor
 from ..math.blend import blend_core_edge_functions
 from ..plasma.equilibrium import load_equilibrium
-from ..plasma.utility import get_subset_name_index
+from ..plasma.utility import get_entry_reference, get_subset_name_index
 
 __all__ = ["load_radiation_emitter"]
 
@@ -135,8 +135,12 @@ def _create_rad_func_ggd(
     grid_ggd: IDSStructure,
     data: NDArray[np.float64],
     grid_subset_id: int,
+    ids_root: IDSStructure,
+    db_args: tuple[Any, ...] | None,
+    db_kwargs: dict[str, Any] | None,
     **interp_kwargs,
 ) -> tuple[AxisymmetricMapper, dict[str, float]]:
+    grid_ggd = _resolve_grid_ggd_reference(grid_ggd, ids_root, db_args, db_kwargs)
     grid, subsets, subset_id = load_grid(grid_ggd, with_subsets=True)
     grid_subset_name, grid_subset_index = get_subset_name_index(subset_id, grid_subset_id)
 
@@ -148,6 +152,49 @@ def _create_rad_func_ggd(
     return rad_func, grid.mesh_extent
 
 
+def _resolve_grid_ggd_reference(
+    grid_ggd: IDSStructure,
+    ids_root: IDSStructure,
+    db_args: tuple[Any, ...] | None,
+    db_kwargs: dict[str, Any] | None,
+) -> IDSStructure:
+    if len(grid_ggd.space):
+        return grid_ggd
+
+    if not len(grid_ggd.path):
+        return grid_ggd
+
+    path = str(grid_ggd.path).strip()
+    if not path:
+        return grid_ggd
+
+    if "#" in path:
+        if db_args is None:
+            raise RuntimeError(
+                "Unable to resolve external grid_ggd.path without DBEntry arguments."
+            )
+        with DBEntry(*db_args, **(db_kwargs or {})) as entry:
+            resolved = load_ids_path_reference(entry, path)
+    else:
+        ids_name = ids_root.metadata.name
+        if path.startswith("/"):
+            resolved = resolve_ids_path_reference(ids_root, f"#{ids_name}{path}")
+        else:
+            resolved = resolve_ids_path_reference(ids_root, path)
+
+    if isinstance(resolved, IDSStructArray):
+        if not len(resolved):
+            raise RuntimeError(f"Resolved grid reference '{path}' points to an empty array.")
+        resolved = resolved[0]
+
+    if not isinstance(resolved, IDSStructure) or not hasattr(resolved, "space"):
+        raise RuntimeError(
+            f"Resolved grid reference '{path}' does not point to a grid_ggd structure."
+        )
+
+    return resolved
+
+
 def load_radiation_emitter(
     *args,
     time: float = 0,
@@ -287,8 +334,8 @@ def load_radiation_emitter(
     # Common variables
     ids = None
     ids2 = None
-    uri: str | None = None
-    uri2: str | None = None
+    entry_reference: str | None = None
+    entry_reference2: str | None = None
     emitter = None
     grid = None
     primitive_name: str | None = None
@@ -306,7 +353,7 @@ def load_radiation_emitter(
                 occurrence=occurrence,
                 time_threshold=time_threshold,
             )
-            uri: str | None = entry.uri
+            entry_reference = get_entry_reference(entry)
     except RuntimeError as err:
         raise RuntimeError("Unable to load radiation IDS.") from err
 
@@ -320,7 +367,7 @@ def load_radiation_emitter(
                     occurrence=occurrence2,
                     time_threshold=time_threshold,
                 )
-                uri2: str | None = entry.uri
+                entry_reference2 = get_entry_reference(entry)
         except RuntimeError as err:
             raise RuntimeError("Unable to load second radiation IDS.") from err
 
@@ -335,6 +382,9 @@ def load_radiation_emitter(
     rad_func = None
     rad_func_core = None
     rad_func_ggd = None
+    ggd_args: tuple[Any, ...] | None = args
+    ggd_kwargs: dict[str, Any] | None = kwargs
+    ggd_ids = ids
 
     if source in {"auto", "values"}:
         # ------------------------------
@@ -371,7 +421,11 @@ def load_radiation_emitter(
                     eq_occurrence = occurrence2
                 if values_ggd is None:
                     values_ggd = values_ggd2
-                uri = f"{uri} + {uri2}"
+                    ggd_args = args2
+                    ggd_kwargs = kwargs2
+                    if ids2 is not None:
+                        ggd_ids = ids2
+                entry_reference = f"{entry_reference} + {entry_reference2}"
 
         if values_core is None and values_ggd is None and source == "values":
             raise RuntimeError(
@@ -407,11 +461,14 @@ def load_radiation_emitter(
             height = zmax - zmin
 
         if values_ggd is not None:
-            grid_ggd = grid_ggd or ids.grid_ggd[0]
+            grid_ggd = grid_ggd or ggd_ids.grid_ggd[0]
             rad_func_ggd, extent = _create_rad_func_ggd(
                 grid_ggd,
                 values_ggd,
                 grid_subset_id,
+                ids_root=ggd_ids,
+                db_args=ggd_args,
+                db_kwargs=ggd_kwargs,
                 interpolator_cache=interpolator_cache,
                 interpolator_cache_dir=interpolator_cache_dir,
             )
@@ -440,16 +497,18 @@ def load_radiation_emitter(
 
         if isinstance(rad_func, Function3D):
             emitter = RadiationFunction(rad_func, step=step)
-            primitive_name = f"RadiationEmitter_{ids.time[0]}s, uri {uri}"
+            primitive_name = f"RadiationEmitter_{ids.time[0]}s, entry {entry_reference}"
 
     # ------------------------------------
     # === Load emissivity coefficients ===
     # ------------------------------------
 
     if emitter is None and source in {"auto", "coefficients"}:
+        grid_ggd = _resolve_grid_ggd_reference(ids.grid_ggd[0], ids, args, kwargs)
+
         # Load GGD Grid
         grid = load_grid(
-            ids.grid_ggd[0],
+            grid_ggd,
             with_subsets=False,
             num_toroidal=num_toroidal,
         )
@@ -489,7 +548,7 @@ def load_radiation_emitter(
                 " no values or coefficients are available."
             )
 
-        constructor = FourierBezierConstructor(ids.grid_ggd[0], coefficients=coeff)
+        constructor = FourierBezierConstructor(grid_ggd, coefficients=coeff)
 
         if phis is None:
             d_phi = 360.0 / grid.num_toroidal
@@ -498,7 +557,7 @@ def load_radiation_emitter(
             phis_array = np.asarray(phis, dtype=np.float64)
 
         emissivity = constructor.average_gaussian_faces_per_toroidal(phis_array).ravel()
-        primitive_name = f"RadiationEmitter_{ids.time[0]}s, uri {uri}"
+        primitive_name = f"RadiationEmitter_{ids.time[0]}s, entry {entry_reference}"
 
         rad_func = grid.interpolator(
             emissivity,
diff --git a/src/cherab/imas/plasma/blend.py b/src/cherab/imas/plasma/blend.py
index 5f131d2..2489e2b 100644
--- a/src/cherab/imas/plasma/blend.py
+++ b/src/cherab/imas/plasma/blend.py
@@ -46,6 +46,7 @@
 from .utility import (
     ZERO_VELOCITY,
     ProfileInterpolator,
+    get_entry_reference,
     get_subset_name_index,
     warn_unsupported_species,
 )
@@ -179,6 +180,7 @@ def load_plasma(
                 occurrence=occurrence_core,
                 time_threshold=time_threshold,
             )
+            core_entry_reference = get_entry_reference(entry_core)
     except RuntimeError:
         return load_edge_plasma(
             *edge_args,
@@ -204,6 +206,7 @@ def load_plasma(
                 occurrence=occurrence_edge,
                 time_threshold=time_threshold,
             )
+            edge_entry_reference = get_entry_reference(entry_edge)
     except RuntimeError:
         return load_core_plasma(
             *args,
@@ -269,7 +272,14 @@ def load_plasma(
 
     # === Edge grid and composition ===
     grid_ggd = grid_ggd or edge_profiles_ids.grid_ggd[0]
-    grid, subsets, subset_id = load_grid(grid_ggd, with_subsets=True)
+    needs_external_grid_reference = (
+        not len(grid_ggd.space) and len(grid_ggd.path) and "#" in str(grid_ggd.path)
+    )
+    if needs_external_grid_reference:
+        with DBEntry(*edge_args, **edge_kwargs) as entry:
+            grid, subsets, subset_id = load_grid(grid_ggd, with_subsets=True, entry=entry)
+    else:
+        grid, subsets, subset_id = load_grid(grid_ggd, with_subsets=True)
 
     try:
         grid_subset_name, grid_subset_index = get_subset_name_index(subset_id, grid_subset_id)
@@ -297,7 +307,7 @@ def load_plasma(
     time_edge = edge_profiles_ids.time[0]
     name = (
         f"IMAS core + edge plasma: core/edge time {time_core}/{time_edge}, "
-        f"uri {entry_core.uri} / {entry_edge.uri}."
+        f"uri {core_entry_reference} / {edge_entry_reference}."
     )
     plasma = Plasma(parent=parent, name=name)
 
diff --git a/src/cherab/imas/plasma/core.py b/src/cherab/imas/plasma/core.py
index 00273a3..c3c7db2 100644
--- a/src/cherab/imas/plasma/core.py
+++ b/src/cherab/imas/plasma/core.py
@@ -39,7 +39,12 @@
 from ..ids.common.grid_radial import get_psi_norm, load_core_grid
 from ..ids.core_profiles import load_core_species
 from .equilibrium import load_equilibrium, load_magnetic_field
-from .utility import ZERO_VELOCITY, ProfileInterpolator, warn_unsupported_species
+from .utility import (
+    ZERO_VELOCITY,
+    ProfileInterpolator,
+    get_entry_reference,
+    warn_unsupported_species,
+)
 
 __all__ = ["load_core_plasma"]
 
@@ -127,6 +132,7 @@ def load_core_plasma(
         core_profiles_ids = get_ids_time_slice(
             entry, "core_profiles", time=time, occurrence=occurrence, time_threshold=time_threshold
         )
+        entry_reference = get_entry_reference(entry)
 
     if not len(core_profiles_ids.profiles_1d):
         raise RuntimeError("The profiles_1d AOS in core_profiles IDS is empty.")
@@ -168,7 +174,7 @@ def load_core_plasma(
     # ----------------------------
     # === Create Plasma object ===
     # ----------------------------
-    name = f"IMAS core plasma: time {core_profiles_ids.time[0]}, uri {entry.uri}."
+    name = f"IMAS core plasma: time {core_profiles_ids.time[0]}, uri {entry_reference}."
     plasma = Plasma(parent=parent, name=name)
     radius_inner, radius_outer = equilibrium.r_range
     zmin, zmax = equilibrium.z_range
diff --git a/src/cherab/imas/plasma/edge.py b/src/cherab/imas/plasma/edge.py
index 67c3b97..c4f8b7c 100644
--- a/src/cherab/imas/plasma/edge.py
+++ b/src/cherab/imas/plasma/edge.py
@@ -46,6 +46,7 @@
 from .utility import (
     ZERO_VELOCITY,
     ProfileInterpolator,
+    get_entry_reference,
     get_subset_name_index,
     warn_unsupported_species,
 )
@@ -130,15 +131,20 @@ def load_edge_plasma(
         edge_profiles_ids = get_ids_time_slice(
             entry, "edge_profiles", time=time, occurrence=occurrence, time_threshold=time_threshold
         )
+        entry_reference = get_entry_reference(entry)
 
-    if not len(edge_profiles_ids.grid_ggd) and grid_ggd is None:
-        raise RuntimeError(
-            "The 'grid_ggd' AOS of the edge_profiles IDS is empty "
-            + "and an alternative grid_ggd structure is not provided."
-        )
+        if not len(edge_profiles_ids.grid_ggd) and grid_ggd is None:
+            raise RuntimeError(
+                "The 'grid_ggd' AOS of the edge_profiles IDS is empty "
+                + "and an alternative grid_ggd structure is not provided."
+            )
+
+        if not len(edge_profiles_ids.ggd):
+            raise RuntimeError("The 'ggd' AOS of the edge_profiles IDS is empty.")
 
-    if not len(edge_profiles_ids.ggd):
-        raise RuntimeError("The 'ggd' AOS of the edge_profiles IDS is empty.")
+        # Resolve path-based grid references while the original entry is still open.
+        grid_ggd_local = grid_ggd or edge_profiles_ids.grid_ggd[0]
+        grid, subsets, subset_id = load_grid(grid_ggd_local, with_subsets=True, entry=entry)
 
     # Load magnetic field data. If not provided, try to load from the equilibrium IDS.
     if b_field is None:
@@ -152,10 +158,6 @@ def load_edge_plasma(
             except RuntimeError:
                 print("Warning! No magnetic field data available in the equilibrium IDS.")
 
-    # Create edge grid
-    grid_ggd = grid_ggd or edge_profiles_ids.grid_ggd[0]
-    grid, subsets, subset_id = load_grid(grid_ggd, with_subsets=True)
-
     try:
         grid_subset_name, grid_subset_index = get_subset_name_index(subset_id, grid_subset_id)
         if not np.array_equal(subsets[grid_subset_name], np.arange(grid.num_cell, dtype=int)):
@@ -178,7 +180,7 @@ def load_edge_plasma(
     # ----------------------------
     # === Create Plasma object ===
     # ----------------------------
-    name = f"IMAS edge plasma: time {edge_profiles_ids.time[0]}, uri {entry.uri}."
+    name = f"IMAS edge plasma: time {edge_profiles_ids.time[0]}, uri {entry_reference}."
     plasma = Plasma(parent=parent, name=name)
 
     # Create plasma geometry
diff --git a/src/cherab/imas/plasma/utility.py b/src/cherab/imas/plasma/utility.py
index ca03af3..93f9f09 100644
--- a/src/cherab/imas/plasma/utility.py
+++ b/src/cherab/imas/plasma/utility.py
@@ -30,6 +30,7 @@
 
 __all__ = [
     "ProfileInterpolator",
+    "get_entry_reference",
     "warn_unsupported_species",
     "get_subset_name_index",
 ]
@@ -38,6 +39,36 @@
 ZERO_VELOCITY = ConstantVector3D(Vector3D(0, 0, 0))
 
 
+def get_entry_reference(entry: object) -> str:
+    """Return a human-readable DBEntry reference.
+
+    Parameters
+    ----------
+    entry
+        DBEntry-like object.
+
+    Returns
+    -------
+    str
+        URI if available, otherwise a fallback string using legacy DBEntry attributes
+        in ``key=value`` format.
+    """
+    uri = getattr(entry, "uri", None)
+    if uri:
+        return str(uri)
+
+    legacy_attr_names = ("backend_id", "db_name", "pulse", "run", "user_name", "data_version")
+    parts = []
+    for name in legacy_attr_names:
+        if hasattr(entry, name):
+            parts.append(f"{name}={getattr(entry, name)!r}")
+
+    if not parts:
+        return "<unknown DBEntry>"
+
+    return ", ".join(parts)
+
+
 @dataclass
 class ProfileInterpolator:
     """Dataclass to hold the interpolators for profiles."""

From 590569ba5b5c447d979cbd01b233969f0f12a594 Mon Sep 17 00:00:00 2001
From: munechika-koyo <munechika.koyo@gmail.com>
Date: Mon, 3 Aug 2026 16:40:29 +0200
Subject: [PATCH 03/36] =?UTF-8?q?=E2=9C=85=20Add=20unit=20tests=20for=20gr?=
 =?UTF-8?q?id=20loading=20and=20entry=20reference=20handling=20in=20plasma?=
 =?UTF-8?q?=20and=20common=20modules?=
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit

---
 tests/ggd/test_load_grid_path_reference.py | 74 ++++++++++++++++++
 tests/ids/common/test_path.py              | 73 ++++++++++++++++++
 tests/plasma/test_utility.py               | 89 ++++++++++++++++++++++
 3 files changed, 236 insertions(+)
 create mode 100644 tests/ggd/test_load_grid_path_reference.py
 create mode 100644 tests/ids/common/test_path.py
 create mode 100644 tests/plasma/test_utility.py

diff --git a/tests/ggd/test_load_grid_path_reference.py b/tests/ggd/test_load_grid_path_reference.py
new file mode 100644
index 0000000..31c3672
--- /dev/null
+++ b/tests/ggd/test_load_grid_path_reference.py
@@ -0,0 +1,74 @@
+from uuid import uuid4
+
+import pytest
+from imas import DBEntry
+from imas.ids_defs import MEMORY_BACKEND
+from imas.ids_factory import IDSFactory
+from imas.ids_toplevel import IDSToplevel
+
+from cherab.imas.ids.common import get_ids_time_slice
+from cherab.imas.ids.common.ggd import load_grid
+
+
+def _build_jintrac_external_grid_reference_entry(
+    path_iter_jintrac: str,
+) -> tuple[DBEntry, IDSToplevel, IDSToplevel]:
+    with DBEntry(path_iter_jintrac, "r") as source_entry:
+        try:
+            edge_profiles = source_entry.get("edge_profiles", autoconvert=False)
+        except Exception as exc:
+            raise RuntimeError(f"JINTRAC dataset does not provide required IDSs: {exc}") from exc
+
+    if not len(edge_profiles.grid_ggd):
+        raise RuntimeError(
+            "JINTRAC dataset does not provide grid_ggd data for this integration test."
+        )
+
+    # Make radiation point to edge_profiles grid and remove its local space.
+    radiation = IDSFactory(edge_profiles._version).new("radiation")
+    radiation.ids_properties.homogeneous_time = edge_profiles.ids_properties.homogeneous_time
+    radiation.time = edge_profiles.time
+    radiation.grid_ggd.resize(1)
+    radiation.grid_ggd[0].space.resize(0)
+    radiation.grid_ggd[0].path = "#edge_profiles/grid_ggd(1)"
+
+    token = uuid4().int
+    entry = DBEntry(
+        MEMORY_BACKEND,
+        f"cherab_grid_ref_{token & 0xFFFF:04x}",
+        1 + token % 2_000_000_000,
+        1 + (token >> 31) % 2_000_000_000,
+        data_version=radiation._version,
+    )
+    entry.create()
+    entry.put(edge_profiles)
+    entry.put(radiation)
+
+    radiation_ids = get_ids_time_slice(entry, "radiation", time=0)
+    edge_ids = get_ids_time_slice(entry, "edge_profiles", time=0)
+    return entry, radiation_ids, edge_ids
+
+
+@pytest.mark.requires_imas_memory_backend
+def test_load_grid_resolves_external_reference_on_jintrac_dataset(path_iter_jintrac: str):
+    entry, radiation_ids, edge_ids = _build_jintrac_external_grid_reference_entry(path_iter_jintrac)
+
+    try:
+        resolved_grid = load_grid(radiation_ids.grid_ggd[0], with_subsets=False, entry=entry)
+        expected_grid = load_grid(edge_ids.grid_ggd[0], with_subsets=False)
+
+        assert type(resolved_grid) is type(expected_grid)
+        assert resolved_grid.num_cell == expected_grid.num_cell
+    finally:
+        entry.close()
+
+
+@pytest.mark.requires_imas_memory_backend
+def test_load_grid_requires_entry_for_external_reference_on_jintrac_dataset(path_iter_jintrac: str):
+    entry, radiation_ids, _ = _build_jintrac_external_grid_reference_entry(path_iter_jintrac)
+
+    try:
+        with pytest.raises(RuntimeError, match="without a DBEntry"):
+            load_grid(radiation_ids.grid_ggd[0], with_subsets=False)
+    finally:
+        entry.close()
diff --git a/tests/ids/common/test_path.py b/tests/ids/common/test_path.py
new file mode 100644
index 0000000..fb489a5
--- /dev/null
+++ b/tests/ids/common/test_path.py
@@ -0,0 +1,73 @@
+import pytest
+from imas import DBEntry
+
+from cherab.imas.ids.common import (
+    IDSPathReference,
+    load_ids_path_reference,
+    parse_ids_path_fragment,
+    resolve_ids_path_reference,
+)
+
+
+@pytest.mark.parametrize(
+    ("reference", "expected"),
+    [
+        (
+            "#core_profiles",
+            IDSPathReference(ids_name="core_profiles"),
+        ),
+        (
+            "#core_profiles:2",
+            IDSPathReference(ids_name="core_profiles", occurrence=2),
+        ),
+        (
+            "imas://example?path=/data#edge_profiles:1/grid_ggd(:)/path",
+            IDSPathReference(
+                ids_name="edge_profiles",
+                occurrence=1,
+                idspath="/grid_ggd(:)/path",
+            ),
+        ),
+        (
+            "#grid_ggd/space(1)/objects_per_dimension(1)",
+            IDSPathReference(
+                ids_name="grid_ggd",
+                idspath="/space(1)/objects_per_dimension(1)",
+            ),
+        ),
+    ],
+)
+def test_parse_ids_path_fragment(reference: str, expected: IDSPathReference) -> None:
+    assert parse_ids_path_fragment(reference) == expected
+
+
+@pytest.mark.parametrize("reference", ["", "#", "#:1", "#core_profiles:", "#core_profiles:abc"])
+def test_parse_ids_path_fragment_rejects_invalid_input(reference: str) -> None:
+    with pytest.raises(ValueError):
+        parse_ids_path_fragment(reference)
+
+
+def test_resolve_ids_path_reference_on_loaded_ids(path_iter_jintrac: str) -> None:
+    from cherab.imas.ids.common import get_ids_time_slice
+
+    with DBEntry(path_iter_jintrac, "r") as entry:
+        ids = get_ids_time_slice(entry, "edge_profiles", time=0)
+
+    reference = IDSPathReference(ids_name=str(ids.metadata.name), idspath="/grid_ggd")
+    resolved = resolve_ids_path_reference(ids, reference)
+
+    assert resolved.metadata.path == ids.grid_ggd.metadata.path
+    assert type(resolved) is type(ids.grid_ggd)
+
+
+def test_load_ids_path_reference_from_entry(path_iter_jintrac: str) -> None:
+    with DBEntry(path_iter_jintrac, "r") as entry:
+        resolved = load_ids_path_reference(entry, "#edge_profiles/grid_ggd")
+
+    from cherab.imas.ids.common import get_ids_time_slice
+
+    with DBEntry(path_iter_jintrac, "r") as entry:
+        ids = get_ids_time_slice(entry, "edge_profiles", time=0)
+
+    assert resolved.metadata.path == ids.grid_ggd.metadata.path
+    assert type(resolved) is type(ids.grid_ggd)
diff --git a/tests/plasma/test_utility.py b/tests/plasma/test_utility.py
new file mode 100644
index 0000000..d2f74cb
--- /dev/null
+++ b/tests/plasma/test_utility.py
@@ -0,0 +1,89 @@
+from contextlib import suppress
+
+import pytest
+from imas import DBEntry
+from imas.ids_defs import MEMORY_BACKEND
+
+from cherab.imas.plasma.utility import get_entry_reference
+
+
+@pytest.mark.parametrize(
+    ("constructor", "entry_kwargs", "expected"),
+    [
+        (
+            "uri",
+            {
+                "uri_builder": lambda _tmp_path, _path_iter_jintrac: (
+                    "imas:memory?path=cherab_test_memory_uri"
+                ),
+                "mode": "w",
+            },
+            "imas:memory?path=cherab_test_memory_uri",
+        ),
+        (
+            "uri",
+            {
+                "uri_builder": lambda _tmp_path, path_iter_jintrac: path_iter_jintrac,
+                "mode": "r",
+            },
+            "__PATH_ITER_JINTRAC__",
+        ),
+        pytest.param(
+            "legacy",
+            {
+                "backend_id": MEMORY_BACKEND,
+                "db_name": "ITER",
+                "pulse": 116100,
+                "run": 1001,
+                "data_version": "3",
+            },
+            f"backend_id={MEMORY_BACKEND!r}, db_name='ITER', pulse=116100, run=1001, "
+            "user_name=None, data_version='3'",
+            marks=pytest.mark.requires_imas_memory_backend,
+            id="memory-legacy-v3",
+        ),
+        pytest.param(
+            "legacy",
+            {
+                "backend_id": MEMORY_BACKEND,
+                "db_name": "ITER",
+                "pulse": 134110,
+                "run": 111,
+                "data_version": "4",
+            },
+            f"backend_id={MEMORY_BACKEND!r}, db_name='ITER', pulse=134110, run=111, "
+            "user_name=None, data_version='4'",
+            marks=pytest.mark.requires_imas_memory_backend,
+            id="memory-legacy-v4",
+        ),
+    ],
+    ids=["memory-uri", "file-uri", None, None],
+)
+def test_get_entry_reference_supports_dbentry_constructor_variants(
+    tmp_path,
+    path_iter_jintrac: str,
+    constructor: str,
+    entry_kwargs: dict,
+    expected: str,
+):
+    """Test get_entry_reference() for both URI and legacy DBEntry constructors."""
+    if constructor == "uri":
+        source_uri = entry_kwargs["uri_builder"](tmp_path, path_iter_jintrac)
+        entry = DBEntry(source_uri, entry_kwargs["mode"])
+    else:
+        entry = DBEntry(
+            entry_kwargs["backend_id"],
+            entry_kwargs["db_name"],
+            entry_kwargs["pulse"],
+            entry_kwargs["run"],
+            data_version=entry_kwargs["data_version"],
+        )
+
+    try:
+        actual = get_entry_reference(entry)
+        expected_value = path_iter_jintrac if expected == "__PATH_ITER_JINTRAC__" else expected
+
+        assert actual == expected_value
+    finally:
+        with suppress(Exception):
+            entry.close()

From d4f57d276160357f2462da9ed7ecb17a08d3c7da Mon Sep 17 00:00:00 2001
From: munechika-koyo <munechika.koyo@gmail.com>
Date: Mon, 3 Aug 2026 16:41:47 +0200
Subject: [PATCH 04/36] =?UTF-8?q?=F0=9F=94=A7=20Remove=20unused=20'uv'=20d?=
 =?UTF-8?q?ependency=20and=20add=20task=20to=20clean=20C=20source=20files?=
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit

---
 pixi.toml | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/pixi.toml b/pixi.toml
index 4dbdded..fd984a5 100644
--- a/pixi.toml
+++ b/pixi.toml
@@ -19,7 +19,6 @@ noarch = false
 python = ["3.10.*", "3.11.*", "3.12.*", "3.13.*", "3.14.*"]
 
 [package.host-dependencies]
-uv = "*"
 python = "*"
 hatch-cython = ">=0.6.0"
 hatch-vcs = "*"
@@ -136,6 +135,7 @@ lint = { cmd = "lefthook run pre-commit --all-files --force", description = "
 # Clean build artifacts
 [feature.tools.tasks]
 clean = { cmd = "find src -type f \\( -name '*.c' -o -name '*.so' -o -name '*.pyd' -o -name '*.dll' \\) -delete", description = "🧹 Clean build artifacts" }
+clean-c = { cmd = "find src -type f -name '*.c' -delete", description = "🧹 Clean C source files" }
 doc-clean = { cmd = [
   "rm",
   "-rf",

From 6e0038e2b32132a698b5ac2e0090e7612fc234f0 Mon Sep 17 00:00:00 2001
From: munechika-koyo <munechika.koyo@gmail.com>
Date: Mon, 3 Aug 2026 17:32:46 +0200
Subject: [PATCH 05/36] =?UTF-8?q?=F0=9F=8E=A8=20Add=20logging=20for=20grid?=
 =?UTF-8?q?=20source=20resolution=20and=20remove=20unused=20lazy=20loading?=
 =?UTF-8?q?=20in=20path=20reference?=
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit

---
 src/cherab/imas/ids/common/ggd/load_grid.py | 2 ++
 src/cherab/imas/ids/common/path.py          | 2 +-
 2 files changed, 3 insertions(+), 1 deletion(-)

diff --git a/src/cherab/imas/ids/common/ggd/load_grid.py b/src/cherab/imas/ids/common/ggd/load_grid.py
index ad27b87..1cdcdf5 100644
--- a/src/cherab/imas/ids/common/ggd/load_grid.py
+++ b/src/cherab/imas/ids/common/ggd/load_grid.py
@@ -174,6 +174,8 @@ def _resolve_grid_source(grid_ggd: IDSStructure, entry: DBEntry | None = None) -
     if not path:
         return grid_ggd
 
+    print(f"Info: resolving GGD grid from grid_ggd.path '{path}'.")
+
     resolved: IDSToplevel | IDSStructure | IDSStructArray
 
     if "#" in path:
diff --git a/src/cherab/imas/ids/common/path.py b/src/cherab/imas/ids/common/path.py
index f0da271..a999bcd 100644
--- a/src/cherab/imas/ids/common/path.py
+++ b/src/cherab/imas/ids/common/path.py
@@ -140,7 +140,7 @@ def load_ids_path_reference(
         The resolved IDS object.
     """
     ref = parse_ids_path_fragment(reference) if isinstance(reference, str) else reference
-    ids = entry.get(ref.ids_name, occurrence=ref.occurrence or 0, autoconvert=False, lazy=True)
+    ids = entry.get(ref.ids_name, occurrence=ref.occurrence or 0, autoconvert=False)
     if not ref.idspath:
         return ids
 

From 031abf624c0a45cdc1e6fb96f1f3c30c84e6a1ba Mon Sep 17 00:00:00 2001
From: munechika-koyo <munechika.koyo@gmail.com>
Date: Mon, 3 Aug 2026 17:33:03 +0200
Subject: [PATCH 06/36] =?UTF-8?q?=F0=9F=93=9D=20Update=20changelog=20for?=
 =?UTF-8?q?=20version=200.5.1:=20add=20IDS=20path=20utility=20helpers,=20u?=
 =?UTF-8?q?nit=20tests,=20and=20enhance=20grid=20loading?=
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit

---
 CHANGELOG.md | 6 +++++-
 1 file changed, 5 insertions(+), 1 deletion(-)

diff --git a/CHANGELOG.md b/CHANGELOG.md
index 13a07cc..18d37a0 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -5,16 +5,20 @@ All notable changes to this project will be documented in this file.
 The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.1.0/),
 and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0.html).
 
-## [0.5.1] - 2026-07-28
+## [0.5.1]
 
 ### Added
 
 - Add synthetic JINTRAC radiation values dataset support and related regression tests
 - Extend radiation emitter loading to support dual emissivity sources with improved validation
+- Add IDS path utility helpers and `get_entry_reference` for resolving entry references
+- Add unit tests for IDS path handling, grid loading via path references, and plasma utility entry-reference workflows
 
 ### Changed
 
 - Improve radiation emitter loading checks for duplicate emissivity values and core-profile grid data
+- Refactor `load_grid` to support explicit `entry` selection and improve referenced-grid source resolution
+- Enhance plasma and emitter loading paths to use entry-reference aware grid resolution
 
 ### Fixed
 

From 1a1c8544d2870305f5520272a330dada92e07100 Mon Sep 17 00:00:00 2001
From: munechika-koyo <munechika.koyo@gmail.com>
Date: Fri, 7 Aug 2026 11:35:57 +0200
Subject: [PATCH 07/36] =?UTF-8?q?=E2=9E=95=20Add=20`rich`=20as=20a=20runti?=
 =?UTF-8?q?me=20dependency?=
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit

---
 CHANGELOG.md   | 1 +
 pixi.toml      | 2 +-
 pyproject.toml | 3 ++-
 3 files changed, 4 insertions(+), 2 deletions(-)

diff --git a/CHANGELOG.md b/CHANGELOG.md
index 18d37a0..de1d9f8 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -19,6 +19,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 - Improve radiation emitter loading checks for duplicate emissivity values and core-profile grid data
 - Refactor `load_grid` to support explicit `entry` selection and improve referenced-grid source resolution
 - Enhance plasma and emitter loading paths to use entry-reference aware grid resolution
+- Promote `rich` from a test-only dependency to a runtime dependency
 
 ### Fixed
 
diff --git a/pixi.toml b/pixi.toml
index fd984a5..d148d73 100644
--- a/pixi.toml
+++ b/pixi.toml
@@ -28,6 +28,7 @@ cherab = "1.5.*"
 
 [package.run-dependencies]
 imas-python = "*"
+rich = "*"
 pooch = "*"
 typing-extensions = ">=4.5"
 
@@ -37,7 +38,6 @@ typing-extensions = ">=4.5"
 [dependencies]
 cherab-imas = { path = "." }
 ipython = "*"
-rich = "*"
 
 # Publication-quality plot
 ultraplot = ">=1.72.0"
diff --git a/pyproject.toml b/pyproject.toml
index 56f240d..958eade 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -34,13 +34,14 @@ classifiers = [
 dependencies = [
   "cherab==1.5.*",
   "imas-python[netcdf]",
+  "rich",
   "pooch",
   "typing-extensions; python_version < '3.12'",
 ]
 dynamic = ["version"]
 
 [project.optional-dependencies]
-test = ["pytest", "pytest-cov", "rich", "plotly"]
+test = ["pytest", "pytest-cov", "plotly"]
 
 [project.urls]
 Homepage = "https://github.com/cherab"

From 522c0e7de51ec15e58164f8474a950c828863f4f Mon Sep 17 00:00:00 2001
From: munechika-koyo <munechika.koyo@gmail.com>
Date: Mon, 10 Aug 2026 11:45:55 +0200
Subject: [PATCH 08/36] =?UTF-8?q?=E2=9C=A8=20Implement=202D=20cell=20geome?=
 =?UTF-8?q?try=20calculation=20for=20triangulated=20polygons?=
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit

---
 src/cherab/imas/math/polygon.pyi |  27 +++++++
 src/cherab/imas/math/polygon.pyx | 127 +++++++++++++++++++++++++++++++
 2 files changed, 154 insertions(+)
 create mode 100644 src/cherab/imas/math/polygon.pyi
 create mode 100644 src/cherab/imas/math/polygon.pyx

diff --git a/src/cherab/imas/math/polygon.pyi b/src/cherab/imas/math/polygon.pyi
new file mode 100644
index 0000000..f268e6c
--- /dev/null
+++ b/src/cherab/imas/math/polygon.pyi
@@ -0,0 +1,27 @@
+from numpy import float64, int32
+from numpy.typing import NDArray
+
+def calculate_2d_cell_geometry(
+    vertices: NDArray[float64],
+    triangles: NDArray[int32],
+    cell_to_triangle: NDArray[int32],
+) -> tuple[NDArray[float64], NDArray[float64]]:
+    """Calculate areas and area centroids of triangulated 2-D cells.
+
+    Each polygonal cell is represented by one or more triangles. The centroid is
+    the area-weighted centroid of those triangles.
+
+    Parameters
+    ----------
+    vertices : (N, 2) ndarray [numpy.float64]
+        Coordinates of all polygon vertices.
+    triangles : (M, 3) ndarray [numpy.int32]
+        Vertex indices of the triangles forming the cells.
+    cell_to_triangle : (K, 2) ndarray [numpy.int32]
+        For each cell, the first triangle index and number of triangles.
+
+    Returns
+    -------
+    tuple[NDArray[float64], NDArray[float64]]
+        Cell centroids with shape ``(K, 2)`` and areas with shape ``(K,)``.
+    """
diff --git a/src/cherab/imas/math/polygon.pyx b/src/cherab/imas/math/polygon.pyx
new file mode 100644
index 0000000..a6b3844
--- /dev/null
+++ b/src/cherab/imas/math/polygon.pyx
@@ -0,0 +1,127 @@
+"""Fast geometry calculations for polygonal meshes."""
+
+import numpy as np
+
+cimport cython
+from cython.parallel import prange
+from numpy cimport import_array, ndarray, int32_t
+
+__all__ = ["calculate_2d_cell_geometry"]
+
+DEF OPENMP_MIN_CELLS = 100000
+
+
+import_array()
+
+
+@cython.boundscheck(False)
+@cython.wraparound(False)
+@cython.cdivision(True)
+cdef void _calculate_cell_geometry(
+    const double[:, ::1] vertices,
+    const int32_t[:, ::1] triangles,
+    const int32_t[:, ::1] cell_to_triangle,
+    Py_ssize_t cell_index,
+    double[:, ::1] cell_centres,
+    double[::1] cell_areas,
+) noexcept nogil:
+    cdef:
+        Py_ssize_t j, triangle_index
+        int32_t vertex0, vertex1, vertex2
+        double x0, y0, x1, y1, x2, y2
+        double triangle_area, area_sum = 0.0
+        double first_moment_x = 0.0, first_moment_y = 0.0
+
+    for j in range(cell_to_triangle[cell_index, 1]):
+        triangle_index = cell_to_triangle[cell_index, 0] + j
+        vertex0 = triangles[triangle_index, 0]
+        vertex1 = triangles[triangle_index, 1]
+        vertex2 = triangles[triangle_index, 2]
+
+        x0, y0 = vertices[vertex0, 0], vertices[vertex0, 1]
+        x1, y1 = vertices[vertex1, 0], vertices[vertex1, 1]
+        x2, y2 = vertices[vertex2, 0], vertices[vertex2, 1]
+
+        triangle_area = 0.5 * abs(
+            (x0 - x2) * (y1 - y2) - (x1 - x2) * (y0 - y2)
+        )
+        area_sum += triangle_area
+        first_moment_x += triangle_area * (x0 + x1 + x2) / 3.0
+        first_moment_y += triangle_area * (y0 + y1 + y2) / 3.0
+
+    cell_areas[cell_index] = area_sum
+    if area_sum > 0.0:
+        cell_centres[cell_index, 0] = first_moment_x / area_sum
+        cell_centres[cell_index, 1] = first_moment_y / area_sum
+    else:
+        cell_centres[cell_index, 0] = 0.0
+        cell_centres[cell_index, 1] = 0.0
+
+
+@cython.boundscheck(False)
+@cython.wraparound(False)
+@cython.initializedcheck(False)
+@cython.cdivision(True)
+cpdef tuple calculate_2d_cell_geometry(
+    const double[:, ::1] vertices,
+    const int32_t[:, ::1] triangles,
+    const int32_t[:, ::1] cell_to_triangle,
+):
+    """Calculate areas and area centroids of triangulated 2-D cells.
+
+    Each polygonal cell is represented by one or more triangles. The centroid is
+    the area-weighted centroid of those triangles.
+
+    Parameters
+    ----------
+    vertices : (N, 2) ndarray [numpy.float64]
+        Coordinates of all polygon vertices.
+    triangles : (M, 3) ndarray [numpy.int32]
+        Vertex indices of the triangles forming the cells.
+    cell_to_triangle : (K, 2) ndarray [numpy.int32]
+        For each cell, the first triangle index and number of triangles.
+
+    Returns
+    -------
+    tuple of ndarray
+        Cell centroids with shape ``(K, 2)`` and areas with shape ``(K,)``.
+    """
+    cdef:
+        Py_ssize_t i
+        ndarray[double, ndim=2] cell_centres
+        ndarray[double, ndim=1] cell_areas
+        double[:, ::1] centres_mv
+        double[::1] areas_mv
+
+    if vertices.shape[1] != 2:
+        raise ValueError("vertices must have a shape of (N, 2).")
+    if triangles.shape[1] != 3:
+        raise ValueError("triangles must have a shape of (M, 3).")
+    if cell_to_triangle.shape[1] != 2:
+        raise ValueError("cell_to_triangle must have a shape of (K, 2).")
+
+    cell_centres = np.empty((cell_to_triangle.shape[0], 2), dtype=np.float64)
+    cell_areas = np.empty(cell_to_triangle.shape[0], dtype=np.float64)
+
+    centres_mv = cell_centres
+    areas_mv = cell_areas
+
+    for i in prange(
+        cell_to_triangle.shape[0],
+        nogil=True,
+        schedule="static",
+        use_threads_if=cell_to_triangle.shape[0] >= OPENMP_MIN_CELLS,
+    ):
+        _calculate_cell_geometry(
+            vertices,
+            triangles,
+            cell_to_triangle,
+            i,
+            centres_mv,
+            areas_mv,
+        )
+
+    if np.any(cell_areas == 0.0):
+        raise ValueError("All cells must have a positive area.")
+
+    return cell_centres, cell_areas

From 562f5cceab80a2e52aa309df44973c00a253b01c Mon Sep 17 00:00:00 2001
From: munechika-koyo <munechika.koyo@gmail.com>
Date: Mon, 10 Aug 2026 11:46:31 +0200
Subject: [PATCH 09/36] =?UTF-8?q?=F0=9F=94=A7=20Fix=20Cython=20compile=20a?=
 =?UTF-8?q?nd=20link=20arguments=20for=20macOS=20support?=
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit

---
 pyproject.toml | 3 ++-
 1 file changed, 2 insertions(+), 1 deletion(-)

diff --git a/pyproject.toml b/pyproject.toml
index 958eade..2a595d9 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -76,7 +76,8 @@ src = "cherab"
 include_numpy = true
 compile_py = false
 parallel = true
-compile_args = ["-O3"]
+compile_args = ["-O3", { arg = "-Xclang=-fopenmp", platforms = ["darwin"] }]
+extra_link_args = [{ arg = "-lomp", platforms = ["darwin"] }]
 env = [
   { env = "LDFLAGS", arg = "-headerpad_max_install_names", platforms = [
     "darwin",

From 6d2a11c9fae6bd95216e107b531759f7f3745e96 Mon Sep 17 00:00:00 2001
From: munechika-koyo <munechika.koyo@gmail.com>
Date: Tue, 11 Aug 2026 08:31:27 +0200
Subject: [PATCH 10/36] =?UTF-8?q?=F0=9F=8E=A8=20Enhance=20grid=20handling?=
 =?UTF-8?q?=20by=20adding=20CellConnectivity,=20CellData=20types=20and=20a?=
 =?UTF-8?q?s=5Fcell=5Fdata=20utility=20function?=
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit

---
 src/cherab/imas/ggd/__init__.py               |  5 +-
 src/cherab/imas/ggd/base_mesh.py              | 77 ++++++++++++++++---
 .../imas/ggd/unstruct_2d_extend_mesh.py       | 41 ++++++----
 src/cherab/imas/ggd/unstruct_3d_mesh.py       | 29 ++++---
 4 files changed, 115 insertions(+), 37 deletions(-)

diff --git a/src/cherab/imas/ggd/__init__.py b/src/cherab/imas/ggd/__init__.py
index fd39c74..d732b40 100644
--- a/src/cherab/imas/ggd/__init__.py
+++ b/src/cherab/imas/ggd/__init__.py
@@ -18,15 +18,18 @@
 """Subpackage for handling general grids (GGD) for IMAS IDSs."""
 
 from . import base_mesh
-from .base_mesh import GGDGrid
+from .base_mesh import CellConnectivity, CellData, GGDGrid, VertexIndices
 from .unstruct_2d_extend_mesh import UnstructGrid2DExtended
 from .unstruct_2d_mesh import UnstructGrid2D
 from .unstruct_3d_mesh import UnstructGrid3D
 
 __all__ = [
     "base_mesh",
+    "CellConnectivity",
+    "CellData",
     "GGDGrid",
     "UnstructGrid2D",
     "UnstructGrid2DExtended",
     "UnstructGrid3D",
+    "VertexIndices",
 ]
diff --git a/src/cherab/imas/ggd/base_mesh.py b/src/cherab/imas/ggd/base_mesh.py
index 8308ff2..0c4ce01 100755
--- a/src/cherab/imas/ggd/base_mesh.py
+++ b/src/cherab/imas/ggd/base_mesh.py
@@ -35,10 +35,27 @@
 from raysect.core.math.function.vector3d.function3d import Function3D as VectorFunction3D
 from raysect.core.math.vector import Vector3D
 
-__all__ = ["GGDGrid", "CellSelection", "InterpolatorCacheMode", "as_index_array"]
+__all__ = [
+    "GGDGrid",
+    "CellConnectivity",
+    "CellData",
+    "CellSelection",
+    "InterpolatorCacheMode",
+    "VertexIndices",
+    "as_cell_data",
+    "as_index_array",
+]
 
 ZEROVECTOR = Vector3D(0, 0, 0)
+
+CellData: TypeAlias = NDArray[np.floating[Any]] | Sequence[float]
+"""One-dimensional scalar data defined on grid cells."""
+VertexIndices: TypeAlias = Sequence[SupportsIndex] | NDArray[np.integer[Any]]
+"""Vertex indices defining one grid cell."""
+CellConnectivity: TypeAlias = NDArray[np.integer[Any]] | Sequence[VertexIndices]
+"""Vertex connectivity defining a sequence of grid cells."""
 CellSelection: TypeAlias = Sequence[SupportsIndex] | NDArray[np.integer[Any]]
+"""One-dimensional array of cell-selection indices."""
 InterpolatorCacheMode: TypeAlias = Literal["none", "memory", "disk"]
 """Cache mode for interpolator templates.
 
@@ -50,6 +67,34 @@
 InterpolatorT = TypeVar("InterpolatorT")
 
 
+def as_cell_data(data: CellData, num_cell: int) -> NDArray[np.float64]:
+    """Return validated scalar cell data as a one-dimensional float array.
+
+    Parameters
+    ----------
+    data
+        Scalar values defined on grid cells.
+    num_cell
+        Expected number of cell values.
+
+    Returns
+    -------
+    NDArray[numpy.float64]
+        Validated one-dimensional cell data.
+
+    Raises
+    ------
+    ValueError
+        If the data is not one-dimensional or does not contain ``num_cell`` values.
+    """
+    data_array = np.asarray_chkfinite(data, dtype=np.float64)
+    if data_array.ndim != 1:
+        raise ValueError("Cell data must be one-dimensional.")
+    if data_array.size != num_cell:
+        raise ValueError(f"Cell data must contain {num_cell} values.")
+    return data_array
+
+
 def as_index_array(indices: CellSelection) -> NDArray[np.intp]:
     """Return cell-selection indices as a NumPy integer array.
 
@@ -162,6 +207,10 @@ def mesh_extent(self) -> dict[str, float]:
     def _interpolator_geometry_hash(self) -> str | None:
         """Return a stable geometry hash based on the grid `vertices` and `cells`.
 
+        ``cells`` may be either a regular NumPy array or a sequence of
+        variable-length index arrays (for example, a mixture of triangular and
+        quadrilateral cells).
+
         Returns
         -------
         str | None
@@ -173,16 +222,26 @@ def _interpolator_geometry_hash(self) -> str | None:
         if vertices is None or cells is None:
             return None
 
-        vertices_array = np.ascontiguousarray(vertices)
-        cells_array = np.ascontiguousarray(cells)
-
         digest = hashlib.blake2b(digest_size=20)
+
+        vertices_array = np.ascontiguousarray(vertices)
         digest.update(str(vertices_array.dtype).encode("ascii"))
         digest.update(np.asarray(vertices_array.shape, dtype=np.int64).tobytes())
         digest.update(vertices_array.tobytes())
-        digest.update(str(cells_array.dtype).encode("ascii"))
-        digest.update(np.asarray(cells_array.shape, dtype=np.int64).tobytes())
-        digest.update(cells_array.tobytes())
+
+        if isinstance(cells, np.ndarray) and cells.dtype != object:
+            cells_array = np.ascontiguousarray(cells)
+            digest.update(str(cells_array.dtype).encode("ascii"))
+            digest.update(np.asarray(cells_array.shape, dtype=np.int64).tobytes())
+            digest.update(cells_array.tobytes())
+        else:
+            digest.update(np.asarray(len(cells), dtype=np.int64).tobytes())
+            for cell in cells:
+                cell_array = np.ascontiguousarray(cell)
+                digest.update(str(cell_array.dtype).encode("ascii"))
+                digest.update(np.asarray(cell_array.shape, dtype=np.int64).tobytes())
+                digest.update(cell_array.tobytes())
+
         return digest.hexdigest()
 
     def _interpolator_cache_key(
@@ -451,7 +510,7 @@ def subset(self, indices: CellSelection, name: str | None = None) -> GGDGrid:
     @abstractmethod
     def interpolator(
         self,
-        grid_data: NDArray[np.float64],
+        grid_data: CellData,
         fill_value: float = 0.0,
         *,
         interpolator_cache: InterpolatorCacheMode = "memory",
@@ -524,7 +583,7 @@ def vector_interpolator(
 
     def plot_mesh(
         self,
-        data: NDArray[np.float64] | None = None,
+        data: CellData | None = None,
         ax: matplotlib.axes.Axes | None = None,
         **grid_styles,
     ) -> matplotlib.axes.Axes:
diff --git a/src/cherab/imas/ggd/unstruct_2d_extend_mesh.py b/src/cherab/imas/ggd/unstruct_2d_extend_mesh.py
index 8572ced..93d2d8f 100644
--- a/src/cherab/imas/ggd/unstruct_2d_extend_mesh.py
+++ b/src/cherab/imas/ggd/unstruct_2d_extend_mesh.py
@@ -38,7 +38,15 @@
 
 from ..math import UnstructGridFunction3D, UnstructGridVectorFunction3D
 from ..math.tetrahedralize import calculate_tetra_volume, cell_to_5tetra
-from .base_mesh import CellSelection, GGDGrid, InterpolatorCacheMode, as_index_array
+from .base_mesh import (
+    CellConnectivity,
+    CellData,
+    CellSelection,
+    GGDGrid,
+    InterpolatorCacheMode,
+    as_cell_data,
+    as_index_array,
+)
 from .unstruct_3d_mesh import UnstructGrid3D
 
 __all__ = ["UnstructGrid2DExtended"]
@@ -77,7 +85,7 @@ class UnstructGrid2DExtended(GGDGrid):
     def __init__(
         self,
         vertices: ArrayLike,
-        cells: ArrayLike,
+        cells: CellConnectivity,
         num_faces: int,
         num_poloidal: int,
         num_toroidal: int,
@@ -86,8 +94,8 @@ def __init__(
     ) -> None:
         vertices = np.array(vertices, dtype=np.float64)
         vertices.setflags(write=False)
-        cells = np.array(cells, dtype=np.int32)
-        cells.setflags(write=False)
+        cells_array: NDArray[np.int32] = np.array(cells, dtype=np.int32)
+        cells_array.setflags(write=False)
 
         if vertices.ndim != 2:
             raise ValueError(
@@ -101,20 +109,20 @@ def __init__(
                 + f"The shape of 'vertices' is {vertices.shape}."
             )
 
-        if cells.ndim != 2:
+        if cells_array.ndim != 2:
             raise ValueError(
                 "Attribute 'cells' must be a 2D array-like. "
-                + f"The number of dimensions in 'cells' is {cells.ndim}."
+                + f"The number of dimensions in 'cells' is {cells_array.ndim}."
             )
 
-        if cells.shape[1] != 8:
+        if cells_array.shape[1] != 8:
             raise ValueError(
                 "Attribute 'cells' must be a (M, 8) array-like. "
-                + f"The shape of 'cells' is {cells.shape}."
+                + f"The shape of 'cells' is {cells_array.shape}."
             )
 
         self._vertices: NDArray[np.float64] = vertices
-        self._cells: NDArray[np.int32] = cells
+        self._cells: NDArray[np.int32] = cells_array
         self._num_faces = num_faces
         self._num_poloidal = num_poloidal
         self._num_toroidal = num_toroidal
@@ -425,7 +433,7 @@ def subset(self, indices: CellSelection, name: str | None = None) -> UnstructGri
     @override
     def interpolator(
         self,
-        grid_data: NDArray[np.float64],
+        grid_data: CellData,
         fill_value: float = 0,
         *,
         interpolator_cache: InterpolatorCacheMode = "memory",
@@ -457,6 +465,7 @@ def interpolator(
         `.UnstructGridFunction3D`
             Interpolator instance.
         """
+        grid_data = as_cell_data(grid_data, self._num_cell)
         return self._build_cached_interpolator(
             interpolator_cls=UnstructGridFunction3D,
             template_builder=lambda: UnstructGridFunction3D(
@@ -570,7 +579,7 @@ def __setstate__(self, state):
         self._initial_setup()
 
     def plot_tetra_mesh(
-        self, data: ArrayLike | None = None, ax: matplotlib.axes.Axes | None = None
+        self, data: CellData | None = None, ax: matplotlib.axes.Axes | None = None
     ) -> None:
         """Plot the tetrahedral mesh grid geometry.
 
@@ -594,7 +603,7 @@ def plot_tetra_mesh(
     @override
     def plot_mesh(
         self,
-        data: ArrayLike | None = None,
+        data: CellData | None = None,
         ax: matplotlib.axes.Axes | None = None,
         **grid_styles,
     ) -> matplotlib.axes.Axes:
@@ -630,7 +639,7 @@ def plot_mesh(
             collection_mesh = PolyCollection(verts, **grid_styles)
         else:
             collection_mesh = PolyCollection(verts)
-            collection_mesh.set_array(data)
+            collection_mesh.set_array(as_cell_data(data, self._num_faces))
         ax.add_collection(collection_mesh)
         ax.set_aspect(1)
         ax.set_xlim(self._mesh_extent["rmin"], self._mesh_extent["rmax"])
@@ -647,7 +656,7 @@ def plot_mesh(
 
     def plot_tri_mesh(
         self,
-        data: ArrayLike,
+        data: CellData,
         ax: matplotlib.axes.Axes | None = None,
         cmap: str = "viridis",
         **kwargs,
@@ -684,9 +693,7 @@ def plot_tri_mesh(
                 f"({self._num_faces=}, {self._num_poloidal=}.)"
             )
 
-        data = np.asarray_chkfinite(data)
-        if data.shape[0] != self._num_faces:
-            raise ValueError("The data array must have the same number of faces as the grid.")
+        data = as_cell_data(data, self._num_faces)
         data = np.repeat(data, 2)
 
         # Create triangulation for the first poloidal plane
diff --git a/src/cherab/imas/ggd/unstruct_3d_mesh.py b/src/cherab/imas/ggd/unstruct_3d_mesh.py
index 08cbf6b..8372477 100644
--- a/src/cherab/imas/ggd/unstruct_3d_mesh.py
+++ b/src/cherab/imas/ggd/unstruct_3d_mesh.py
@@ -33,7 +33,15 @@
 
 from ..math import UnstructGridFunction3D, UnstructGridVectorFunction3D
 from ..math.tetrahedralize import calculate_tetra_volume, cell_to_5tetra
-from .base_mesh import CellSelection, GGDGrid, InterpolatorCacheMode, as_index_array
+from .base_mesh import (
+    CellConnectivity,
+    CellData,
+    CellSelection,
+    GGDGrid,
+    InterpolatorCacheMode,
+    as_cell_data,
+    as_index_array,
+)
 
 __all__ = ["UnstructGrid3D"]
 
@@ -61,13 +69,13 @@ class UnstructGrid3D(GGDGrid):
     def __init__(
         self,
         vertices: ArrayLike,
-        cells: ArrayLike,
+        cells: CellConnectivity,
         name: str = "Cells",
     ) -> None:
         vertices = np.array(vertices, dtype=np.float64)
         vertices.setflags(write=False)
-        cells = np.array(cells, dtype=np.int32)
-        cells.setflags(write=False)
+        cells_array: NDArray[np.int32] = np.array(cells, dtype=np.int32)
+        cells_array.setflags(write=False)
 
         if vertices.ndim != 2:
             raise ValueError(
@@ -81,20 +89,20 @@ def __init__(
                 + f"The shape of 'vertices' is {vertices.shape}."
             )
 
-        if cells.ndim != 2:
+        if cells_array.ndim != 2:
             raise ValueError(
                 "Attribute 'cells' must be a 2D array-like. "
-                + f"The number of dimensions in 'cells' is {cells.ndim}."
+                + f"The number of dimensions in 'cells' is {cells_array.ndim}."
             )
 
-        if cells.shape[1] != 8:
+        if cells_array.shape[1] != 8:
             raise ValueError(
                 "Attribute 'cells' must be a (M, 8) array-like. "
-                + f"The shape of 'cells' is {cells.shape}."
+                + f"The shape of 'cells' is {cells_array.shape}."
             )
 
         self._vertices: NDArray[np.float64] = vertices
-        self._cells: NDArray[np.int32] = cells
+        self._cells: NDArray[np.int32] = cells_array
 
         super().__init__(name=name, dimension=3, coordinate_system="cartesian")
 
@@ -252,7 +260,7 @@ def subset(self, indices: CellSelection, name: str | None = None) -> UnstructGri
     @override
     def interpolator(
         self,
-        grid_data: NDArray[np.float64],
+        grid_data: CellData,
         fill_value: float = 0,
         *,
         interpolator_cache: InterpolatorCacheMode = "memory",
@@ -284,6 +292,7 @@ def interpolator(
         `.UnstructGridFunction3D`
             Interpolator instance.
         """
+        grid_data = as_cell_data(grid_data, self._num_cell)
         return self._build_cached_interpolator(
             interpolator_cls=UnstructGridFunction3D,
             template_builder=lambda: UnstructGridFunction3D(

From 2d8106285974a75308d1a5d6135bc8bee6369566 Mon Sep 17 00:00:00 2001
From: munechika-koyo <munechika.koyo@gmail.com>
Date: Tue, 11 Aug 2026 08:41:57 +0200
Subject: [PATCH 11/36] =?UTF-8?q?=E2=9C=A8=20Refactor=20UnstructGrid2D:=20?=
 =?UTF-8?q?Enhance=20valid=20data=20handling=20and=20integrate=20cell=20ge?=
 =?UTF-8?q?ometry=20calculations?=
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit

---
 src/cherab/imas/ggd/unstruct_2d_mesh.py | 199 ++++++++++++++++--------
 1 file changed, 137 insertions(+), 62 deletions(-)

diff --git a/src/cherab/imas/ggd/unstruct_2d_mesh.py b/src/cherab/imas/ggd/unstruct_2d_mesh.py
index 91abb86..3ac278c 100755
--- a/src/cherab/imas/ggd/unstruct_2d_mesh.py
+++ b/src/cherab/imas/ggd/unstruct_2d_mesh.py
@@ -19,8 +19,8 @@
 
 from __future__ import annotations
 
-import hashlib
 import sys
+from collections.abc import Sequence
 from pathlib import Path
 from typing import Any, Literal, cast
 
@@ -38,13 +38,56 @@
 from raysect.core.math.vector import Vector3D
 
 from ..math import UnstructGridFunction2D, UnstructGridVectorFunction2D
-from .base_mesh import CellSelection, GGDGrid, InterpolatorCacheMode, as_index_array
+from ..math.polygon import calculate_2d_cell_geometry
+from .base_mesh import (
+    CellConnectivity,
+    CellData,
+    CellSelection,
+    GGDGrid,
+    InterpolatorCacheMode,
+    as_index_array,
+)
 
 __all__ = ["UnstructGrid2D"]
 
 ZERO_VECTOR = Vector3D(0, 0, 0)
 
 
+def _as_cell_data(data: CellData, valid_data_mask: NDArray[np.bool_]) -> NDArray[np.float64]:
+    """Return validated scalar cell data as a one-dimensional float array.
+
+    Parameters
+    ----------
+    data
+        Scalar values defined on grid cells.
+    valid_data_mask
+        Boolean array indicating valid cell data.
+
+    Returns
+    -------
+    NDArray[numpy.float64]
+        Validated one-dimensional cell data.
+
+    Raises
+    ------
+    ValueError
+        If the data is not one-dimensional or does not match the valid data mask.
+    """
+    data_array = np.asarray_chkfinite(data, dtype=np.float64)
+    if data_array.ndim != 1:
+        raise ValueError("Cell data must be one-dimensional.")
+    num_valid = int(np.count_nonzero(valid_data_mask))
+    if data_array.size == num_valid:
+        return data_array
+    if data_array.size == valid_data_mask.size:
+        return data_array[valid_data_mask]
+    raise ValueError(
+        "Cell data must contain either the number of valid cells or the number of "
+        f"source entries. Data size: {data_array.size}, valid cells: {num_valid}, "
+        f"source entries: {valid_data_mask.size}."
+    )
+
+
 class UnstructGrid2D(GGDGrid):
     """Unstructured 2D grid object.
 
@@ -57,9 +100,14 @@ class UnstructGrid2D(GGDGrid):
     vertices
         Array-like of shape ``(N, 2)`` containing coordinates of the polygon vertices.
     cells
-        List of ``(N,)``-shaped arrays containing the vertex indices in clockwise or
-        counterclockwise order for each polygonal cell in the list (the starting vertex must not be
-        included twice).
+        An ``(N, 4)`` integer array or a list/tuple of 1-D integer index arrays containing
+        the vertex indices in clockwise or counterclockwise order for each polygonal cell
+        (the starting vertex must not be included twice).
+    valid_data_mask
+        Boolean mask over the source face/data entries. Its number of ``True``
+        values must equal the number of cells retained by this grid. Data passed
+        to plotting/interpolation may therefore be either source-sized or already
+        compacted to the valid cells.
     name
         Name of the grid, by default ``'Cells'``.
     coordinate_system
@@ -69,11 +117,12 @@ class UnstructGrid2D(GGDGrid):
     def __init__(
         self,
         vertices: ArrayLike,
-        cells: list[ArrayLike],
+        cells: CellConnectivity,
+        valid_data_mask: NDArray[np.bool_] | Sequence[bool] | None = None,
         name: str = "Cells",
         coordinate_system: Literal["cylindrical", "cartesian"] = "cylindrical",
     ) -> None:
-        vertices = np.asarray_chkfinite(vertices, dtype=np.float64)
+        vertices = np.ascontiguousarray(np.asarray_chkfinite(vertices, dtype=np.float64))
         vertices.setflags(write=False)
 
         if vertices.ndim != 2:
@@ -101,6 +150,20 @@ def __init__(
         self._vertices: NDArray[np.float64] = vertices
         self._cells: tuple[NDArray[np.intp], ...] = tuple(normalized_cells)
 
+        if valid_data_mask is None:
+            valid_data_mask = np.ones(len(self._cells), dtype=np.bool_)
+        else:
+            valid_data_mask = np.asarray(valid_data_mask, dtype=np.bool_)
+        if valid_data_mask.ndim != 1:
+            raise ValueError("valid_data_mask must be one-dimensional.")
+        if (n := np.count_nonzero(valid_data_mask)) != len(self._cells):
+            raise ValueError(
+                f"The number of valid data mask entries ({n})"
+                f" must match the number of cells ({len(self._cells)})."
+            )
+        self._valid_data_mask = np.array(valid_data_mask, dtype=np.bool_, copy=True)
+        self._valid_data_mask.setflags(write=False)
+
         super().__init__(name, 2, coordinate_system)
 
     @override
@@ -138,7 +201,7 @@ def _initial_setup(self) -> None:
         for i, cell in enumerate(self._cells):
             ntri = len(cell) - 2
             if ntri == 1:
-                self._triangles[i] = cell
+                self._triangles[itri] = cell
             else:
                 vert = self._vertices[cell]
                 tri = triangulate2d(cast(Any, vert))
@@ -151,27 +214,16 @@ def _initial_setup(self) -> None:
         self._cell_to_triangle_map.setflags(write=False)
         self._triangle_to_cell_map.setflags(write=False)
 
-        # Calculate cell area and centroid
-        self._cell_centre = np.empty((len(self._cells), 2), dtype=np.float64)
-        self._cell_area = np.empty(len(self._cells), dtype=np.float64)
-
-        vx = x[self._triangles]
-        vy = y[self._triangles]
-        area = 0.5 * np.abs(
-            (vx[:, 0] - vx[:, 2]) * (vy[:, 1] - vy[:, 2])
-            - (vx[:, 1] - vx[:, 2]) * (vy[:, 0] - vy[:, 2])
+        # Calculate cell areas and area centroids in Cython.
+        self._cell_centre, self._cell_area = calculate_2d_cell_geometry(
+            self._vertices, self._triangles, self._cell_to_triangle_map
         )
-
-        for i, cell in enumerate(self._cells):
-            self._cell_centre[i] = self._vertices[cell].mean(0)
-            i_start, ntri = self._cell_to_triangle_map[i]
-            self._cell_area[i] = area[i_start : i_start + ntri].sum()
-
         self._cell_centre.setflags(write=False)
         self._cell_area.setflags(write=False)
 
         if self._coordinate_system == "cylindrical":
-            self._cell_volume = 0.5 * np.pi * self._cell_centre[:, 0] * self._cell_area
+            self._cell_volume = np.multiply(self._cell_centre[:, 0], self._cell_area)
+            np.multiply(self._cell_volume, 2.0 * np.pi, out=self._cell_volume)
             self._cell_volume.setflags(write=False)
 
     @property
@@ -184,34 +236,6 @@ def cells(self) -> tuple[NDArray[np.intp], ...]:
         """List of ``K`` polygonal cells as 1-D integer index arrays."""
         return self._cells
 
-    @override
-    def _interpolator_geometry_hash(self) -> str | None:
-        """Return a stable geometry hash for polygonal 2-D grids.
-
-        This override handles ragged polygon connectivity (`tuple` of variable-length
-        arrays), which cannot be hashed robustly via a single contiguous array.
-
-        Returns
-        -------
-        str | None
-            Stable digest string for cache keys.
-        """
-        digest = hashlib.blake2b(digest_size=20)
-
-        vertices_array = np.ascontiguousarray(self._vertices)
-        digest.update(str(vertices_array.dtype).encode("ascii"))
-        digest.update(np.asarray(vertices_array.shape, dtype=np.int64).tobytes())
-        digest.update(vertices_array.tobytes())
-
-        digest.update(np.asarray(len(self._cells), dtype=np.int64).tobytes())
-        for cell in self._cells:
-            cell_array = np.ascontiguousarray(cell, dtype=np.intp)
-            digest.update(str(cell_array.dtype).encode("ascii"))
-            digest.update(np.asarray(cell_array.shape, dtype=np.int64).tobytes())
-            digest.update(cell_array.tobytes())
-
-        return digest.hexdigest()
-
     @property
     def triangles(self) -> NDArray[np.int32]:
         """Mesh triangles as ``(M, 3)`` array."""
@@ -234,14 +258,27 @@ def cell_to_triangle_map(self) -> NDArray[np.int32]:
         """
         return self._cell_to_triangle_map
 
+    @property
+    def valid_data_mask(self) -> NDArray[np.bool_]:
+        """Boolean mask over source data entries retained by this grid."""
+        return self._valid_data_mask
+
     @override
-    def subset(self, indices: CellSelection, name: str | None = None) -> UnstructGrid2D:
+    def subset(
+        self,
+        indices: CellSelection,
+        name: str | None = None,
+        *,
+        valid_data_mask: NDArray[np.bool_] | Sequence[bool] | None = None,
+    ) -> UnstructGrid2D:
         """Create a subset UnstructGrid2D from this instance.
 
         Parameters
         ----------
         indices
             Indices of the cells of the original grid in the subset.
+        valid_data_mask
+            Boolean array indicating which cells in the subset have valid data.
         name
             Name of the grid subset. Default is ``instance.name + " subset"``.
 
@@ -249,9 +286,34 @@ def subset(self, indices: CellSelection, name: str | None = None) -> UnstructGri
         -------
         `.UnstructGrid2D`
             Subset instance.
+
+        Raises
+        ------
+        ValueError
+            If the validity mask is not one-dimensional or does not select exactly
+            one valid entry per subset cell.
         """
+        # ``load_unstruct_grid_2d(..., with_subsets=True)`` returns the index
+        # array and its source validity mask together. Accept that pair directly
+        # for convenience while retaining the normal ``indices`` API.
+        if valid_data_mask is None and isinstance(indices, tuple) and len(indices) == 2:
+            candidate_mask = np.asarray(indices[1])
+            if candidate_mask.ndim == 1 and candidate_mask.dtype == np.bool_:
+                indices, valid_data_mask = cast(Any, indices[0]), candidate_mask
+
         index_array = as_index_array(indices)
 
+        if valid_data_mask is None:
+            valid_data_mask = np.ones(index_array.size, dtype=np.bool_)
+        else:
+            valid_data_mask = np.asarray(valid_data_mask, dtype=np.bool_)
+        if valid_data_mask.ndim != 1:
+            raise ValueError("valid_data_mask must be one-dimensional.")
+        if np.count_nonzero(valid_data_mask) != index_array.size:
+            raise ValueError(
+                "The number of valid data entries must match the number of subset cells."
+            )
+
         grid = UnstructGrid2D.__new__(UnstructGrid2D)
 
         grid._name = name or self.name + " subset"
@@ -259,6 +321,8 @@ def subset(self, indices: CellSelection, name: str | None = None) -> UnstructGri
         grid._dimension = self._dimension
         grid._scalar_interpolator = None
         grid._vector_interpolator = None
+        grid._valid_data_mask = np.array(valid_data_mask, dtype=np.bool_, copy=True)
+        grid._valid_data_mask.setflags(write=False)
 
         index_list = [int(i) for i in index_array]
         cells_original: tuple[NDArray[np.intp], ...] = tuple(
@@ -282,13 +346,16 @@ def subset(self, indices: CellSelection, name: str | None = None) -> UnstructGri
             i_start += num_vertices
         grid._cells = tuple(cells)
         grid._num_cell = len(grid._cells)
-        ntri_total = i_start - 2 * len(cells_original)
+        ntri_total = sum(len(cell) - 2 for cell in cells)
 
         # cell area and centres of this subset
         grid._cell_area = np.array(self.cell_area[index_array])
         grid._cell_area.setflags(write=False)
         grid._cell_centre = np.array(self.cell_centre[index_array])
         grid._cell_centre.setflags(write=False)
+        if self._coordinate_system == "cylindrical":
+            grid._cell_volume = np.array(self.cell_volume[index_array])
+            grid._cell_volume.setflags(write=False)
 
         # mesh extent of this subset
         xmin, ymin = grid._vertices.min(0)
@@ -318,7 +385,7 @@ def subset(self, indices: CellSelection, name: str | None = None) -> UnstructGri
         for i, cell in enumerate(cells):
             ntri = len(cell) - 2
             if ntri == 1:
-                grid._triangles[i] = cell
+                grid._triangles[itri] = cell
             else:
                 c2t = c2t_map[i]
                 tri = self.triangles[c2t[0] : c2t[0] + c2t[1]]
@@ -336,7 +403,7 @@ def subset(self, indices: CellSelection, name: str | None = None) -> UnstructGri
     @override
     def interpolator(
         self,
-        grid_data: NDArray[np.float64],
+        grid_data: CellData,
         fill_value: float = 0,
         *,
         interpolator_cache: InterpolatorCacheMode = "memory",
@@ -368,6 +435,7 @@ def interpolator(
         `.UnstructGridFunction2D`
             Interpolator instance.
         """
+        grid_data = _as_cell_data(grid_data, self._valid_data_mask)
         return self._build_cached_interpolator(
             interpolator_cls=UnstructGridFunction2D,
             template_builder=lambda: UnstructGridFunction2D(
@@ -455,6 +523,7 @@ def __getstate__(self):
             "coordinate_system": self._coordinate_system,
             "vertices": self._vertices,
             "cells": self._cells,
+            "valid_data_mask": self._valid_data_mask,
         }
         return state
 
@@ -466,12 +535,17 @@ def __setstate__(self, state):
         self._vertices = state["vertices"]
         self._vertices.setflags(write=False)
         self._cells = tuple(np.asarray(cell, dtype=np.intp) for cell in state["cells"])
+        self._valid_data_mask = np.asarray(
+            state.get("valid_data_mask", np.ones(len(self._cells), dtype=np.bool_)),
+            dtype=np.bool_,
+        )
+        self._valid_data_mask.setflags(write=False)
 
         self._initial_setup()
 
     def plot_triangle_mesh(
         self,
-        data: ArrayLike | None = None,
+        data: CellData | None = None,
         ax: matplotlib.axes.Axes | None = None,
         **grid_styles,
     ) -> matplotlib.axes.Axes:
@@ -501,11 +575,12 @@ def plot_triangle_mesh(
         grid_styles.setdefault("linewidth", 0.25)
 
         verts = self._vertices[self._triangles]
+        polygons = cast(Sequence[ArrayLike], verts)
         if data is None:
-            collection_mesh = PolyCollection([verts], **grid_styles)
+            collection_mesh = PolyCollection(polygons, **grid_styles)
         else:
-            data_array = np.asarray(data)
-            collection_mesh = PolyCollection([verts])
+            data_array = _as_cell_data(data, self._valid_data_mask)
+            collection_mesh = PolyCollection(polygons)
             collection_mesh.set_array(data_array[self._triangle_to_cell_map])
         ax.add_collection(collection_mesh)
         ax.set_aspect(1)
@@ -524,7 +599,7 @@ def plot_triangle_mesh(
     @override
     def plot_mesh(
         self,
-        data: ArrayLike | None = None,
+        data: CellData | None = None,
         ax: matplotlib.axes.Axes | None = None,
         **grid_styles,
     ) -> matplotlib.axes.Axes:
@@ -560,7 +635,7 @@ def plot_mesh(
             collection_mesh = PolyCollection(verts, **grid_styles)
         else:
             collection_mesh = PolyCollection(verts)
-            collection_mesh.set_array(data)
+            collection_mesh.set_array(_as_cell_data(data, self._valid_data_mask))
         ax.add_collection(collection_mesh)
         ax.set_aspect(1)
         ax.set_xlim(self._mesh_extent["xmin"], self._mesh_extent["xmax"])

From 1da7cf6316b6270f6f38d6924d173cedca452823 Mon Sep 17 00:00:00 2001
From: munechika-koyo <munechika.koyo@gmail.com>
Date: Tue, 11 Aug 2026 08:43:13 +0200
Subject: [PATCH 12/36] =?UTF-8?q?=F0=9F=8E=A8=20Update=20load=5Fgrid=20fun?=
 =?UTF-8?q?ction=20signature=20to=20include=20boolean=20array=20in=20retur?=
 =?UTF-8?q?n=20type?=
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit

---
 src/cherab/imas/ids/common/ggd/load_grid.py | 5 +++--
 1 file changed, 3 insertions(+), 2 deletions(-)

diff --git a/src/cherab/imas/ids/common/ggd/load_grid.py b/src/cherab/imas/ids/common/ggd/load_grid.py
index 1cdcdf5..88ea774 100644
--- a/src/cherab/imas/ids/common/ggd/load_grid.py
+++ b/src/cherab/imas/ids/common/ggd/load_grid.py
@@ -21,6 +21,7 @@
 
 from typing import Literal, cast, overload
 
+import numpy as np
 from numpy import int32
 from numpy.typing import NDArray
 
@@ -55,7 +56,7 @@ def load_grid(
     num_toroidal: int | None = None,
     *,
     entry: DBEntry | None = None,
-) -> tuple[UnstructGrid2D, dict[str, NDArray[int32]], dict[str, int]]: ...
+) -> tuple[UnstructGrid2D, dict[str, tuple[NDArray[int32], NDArray[np.bool_]]], dict[str, int]]: ...
 
 
 def load_grid(
@@ -66,7 +67,7 @@ def load_grid(
     entry: DBEntry | None = None,
 ) -> (
     UnstructGrid2D
-    | tuple[UnstructGrid2D, dict[str, NDArray[int32]], dict[str, int]]
+    | tuple[UnstructGrid2D, dict[str, tuple[NDArray[int32], NDArray[np.bool_]]], dict[str, int]]
     | UnstructGrid2DExtended
 ):
     """Load grid from the ``grid_ggd`` structure.

From d4c0475bc1eca68f5f9d5baa9ba5b1e3e7fd0992 Mon Sep 17 00:00:00 2001
From: munechika-koyo <munechika.koyo@gmail.com>
Date: Tue, 11 Aug 2026 09:01:30 +0200
Subject: [PATCH 13/36] =?UTF-8?q?=E2=9C=A8=20Enhance=20load=5Funstruct=5Fg?=
 =?UTF-8?q?rid=5F2d:=20Update=20return=20type=20to=20include=20validity=20?=
 =?UTF-8?q?information=20for=20grid=20subsets?=
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit

---
 .../imas/ids/common/ggd/load_unstruct_2d.py   | 92 +++++++++++++------
 1 file changed, 66 insertions(+), 26 deletions(-)

diff --git a/src/cherab/imas/ids/common/ggd/load_unstruct_2d.py b/src/cherab/imas/ids/common/ggd/load_unstruct_2d.py
index a2ff725..a010d4c 100644
--- a/src/cherab/imas/ids/common/ggd/load_unstruct_2d.py
+++ b/src/cherab/imas/ids/common/ggd/load_unstruct_2d.py
@@ -18,18 +18,24 @@
 """Module for loading unstructured 2D grids from IMAS grid_ggd IDS structure."""
 
 from enum import IntEnum
-from typing import Literal, overload
+from typing import Final, Literal, cast, overload
 
 import numpy as np
 from numpy.typing import NDArray
 
-from imas.ids_defs import EMPTY_INT
-from imas.ids_structure import IDSStructure
+from imas.ids_structure import IDSStructArray, IDSStructure
 
 from ....ggd.unstruct_2d_mesh import UnstructGrid2D
 
 __all__ = ["load_unstruct_grid_2d"]
 
+_KNOWN_2D_SUBSET_IDS: Final[frozenset[int]] = frozenset({5, 22, 23, 24, 25, 38, 39, 40})
+"""Set of known 2D grid subset indices in IMAS grid_ggd structure.
+
+The detailed description of each subset can be found in the IMAS data dictionary:
+https://imas-data-dictionary.readthedocs.io/en/latest/generated/identifier/ggd_subset_identifier.html
+"""
+
 
 class DIMENSION(IntEnum):
     """Enumeration for grid dimensions."""
@@ -51,12 +57,17 @@ def load_unstruct_grid_2d(
     space_index: int = 0,
     *,
     with_subsets: Literal[True],
-) -> tuple[UnstructGrid2D, dict[str, NDArray[np.int32]], dict[str, int]]: ...
+) -> tuple[
+    UnstructGrid2D, dict[str, tuple[NDArray[np.int32], NDArray[np.bool_]]], dict[str, int]
+]: ...
 
 
 def load_unstruct_grid_2d(
     grid_ggd: IDSStructure, space_index: int = 0, with_subsets: bool = False
-) -> UnstructGrid2D | tuple[UnstructGrid2D, dict[str, NDArray[np.int32]], dict[str, int]]:
+) -> (
+    UnstructGrid2D
+    | tuple[UnstructGrid2D, dict[str, tuple[NDArray[np.int32], NDArray[np.bool_]]], dict[str, int]]
+):
     """Load unstructured 2D grid from the grid_ggd structure.
 
     Parameters
@@ -72,9 +83,10 @@ def load_unstruct_grid_2d(
     -------
     grid : `.UnstructGrid2D`
         Unstructured 2D grid object.
-    subsets : `dict[str, NDArray[numpy.int32]]`
-        Dictionary with grid subsets for each subset name containing the indices of the cells from
-        that subset. Note that 'Cells' subset is included only if cell indices are specified.
+    subsets : `dict[str, tuple[NDArray[np.int32], NDArray[np.bool_]]]`
+        Dictionary with grid subsets for each subset name containing a tuple with the indices of the
+        cells from that subset and a boolean array indicating the validity of each index.
+        Note that 'Cells' subset is included only if cell indices are specified.
     subset_id : `dict[str, int]`
         Dictionary with grid subset indices.
 
@@ -98,15 +110,32 @@ def load_unstruct_grid_2d(
         vertices[i] = space.objects_per_dimension[DIMENSION.VERTEX].object[i].geometry[:2]
 
     # Reading polygonal cells
+    faces = cast(IDSStructArray, space.objects_per_dimension[DIMENSION.FACE].object)
+    num_faces = len(faces)
     cells = []
+    # ``cells`` is compact (invalid faces are omitted), while GGD subset
+    # references are expressed in the original face numbering.
+    face_to_cell = np.full(num_faces, -1, dtype=np.int32)
+    valid_face = np.ones(num_faces, dtype=bool)
     winding_ok = True
-    for object in space.objects_per_dimension[DIMENSION.FACE].object:
-        # getting cell from nodes
-        cell = np.asarray_chkfinite(object.nodes, dtype=np.int32) - 1  # Fortran to C indexing
+
+    for i_face in range(num_faces):
+        face = faces[i_face]
+
+        if not face.has_value or not face.nodes.has_value:
+            valid_face[i_face] = False
+            continue
+        # Convert every face from Fortran to C indexing. Triangular faces are already ordered.
+        # Only polygons need their winding reconstructed below.
+        cell = np.asarray_chkfinite(face.nodes, dtype=np.int32) - 1
+        if cell.size < 3:
+            valid_face[i_face] = False
+            continue
+
         if cell.size > 3:
             # trying to get the nodes in winding order by parsing the edges
-            edge_dict = {}
-            for boundary in object.boundary:
+            edge_dict: dict[int, list[int]] = {}
+            for boundary in face.boundary:
                 n1, n2 = (
                     space.objects_per_dimension[DIMENSION.EDGE].object[boundary.index - 1].nodes - 1
                 )  # Fortran to C indexing
@@ -127,6 +156,7 @@ def load_unstruct_grid_2d(
                     edge_dict[n2][1] = n1
                 else:
                     edge_dict[n2] = [n1, -1]
+
             if len(edge_dict) == cell.size:  # success, getting the cell nodes in winding order
                 cell1 = np.empty(len(edge_dict), dtype=np.int32)
                 cell1[0] = cell[0]
@@ -139,36 +169,46 @@ def load_unstruct_grid_2d(
             else:
                 winding_ok = False
 
-            cells.append(cell)
+        face_to_cell[i_face] = len(cells)
+        cells.append(cell)
 
     if not winding_ok:
         print("Warning! Unable to verify that the cell nodes are in the winding order.")
 
-    grid = UnstructGrid2D(vertices, cells, name=grid_name)
+    grid = UnstructGrid2D(vertices, cells, valid_face, name=grid_name)
 
     if not with_subsets:
         return grid
 
     # Reading grid subsets (2D only)
-    CELL_SUBSET_IDS = {5, 22, 23, 24, 25, 38, 39, 40}
-    subsets = {}
-    subset_id = {}
+    subsets: dict[str, tuple[NDArray[np.int32], NDArray[np.bool_]]] = {}
+    subset_id: dict[str, int] = {}
     for subset in grid_ggd.grid_subset:
-        subset_index = subset.identifier.index.value
-        dimension_is_2d = subset.dimension == DIMENSION.FACE + 1  # C to Fortran indexing
-        known_subset_id = subset.dimension != EMPTY_INT and subset_index in CELL_SUBSET_IDS
+        subset_index: int = subset.identifier.index.value
+        dimension_is_2d: bool = subset.dimension == DIMENSION.FACE + 1  # C to Fortran indexing
+        known_subset_id: bool = subset_index in _KNOWN_2D_SUBSET_IDS
         if (dimension_is_2d or known_subset_id) and len(subset.element):
             name = str(subset.identifier.name)
-            indices = np.empty(len(subset.element), dtype=np.int32)
-            for i, element in enumerate(subset.element):
+            num_elm = len(subset.element)
+            indices = np.empty(num_elm, dtype=np.int32)
+            valid_subset = np.ones_like(indices, dtype=bool)
+            for i_elm, element in enumerate(subset.element):
                 if len(element.object) > 1:
                     print(
                         f"Warning! Skipping grid subset {name}, "
                         + "because it includes cells not present in the original grid."
                     )
                     break
-                indices[i] = element.object[0].index.value
-            subsets[name] = indices - 1  # Fortran to C indexing
-            subset_id[name] = subset_index
+                face_index = element.object[0].index.value - 1  # Fortran to C indexing
+                if face_index < 0 or face_index >= num_faces:
+                    valid_subset[i_elm] = False
+                    indices[i_elm] = -1
+                    continue
+                indices[i_elm] = face_to_cell[face_index]
+                if indices[i_elm] < 0:
+                    valid_subset[i_elm] = False
+            else:
+                subsets[name] = (indices[valid_subset], valid_subset)
+                subset_id[name] = subset_index
 
     return grid, subsets, subset_id

From 4cef96ea8734229346235916bdf5af6a11e557b5 Mon Sep 17 00:00:00 2001
From: munechika-koyo <munechika.koyo@gmail.com>
Date: Tue, 11 Aug 2026 09:02:34 +0200
Subject: [PATCH 14/36] =?UTF-8?q?=F0=9F=94=A7=20Fix=20num=5Fcells=20calcul?=
 =?UTF-8?q?ation=20to=20correctly=20reference=20the=20first=20element=20of?=
 =?UTF-8?q?=20the=20subset?=
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit

---
 src/cherab/imas/datasets/_builtin.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/cherab/imas/datasets/_builtin.py b/src/cherab/imas/datasets/_builtin.py
index 5f56daa..f233941 100644
--- a/src/cherab/imas/datasets/_builtin.py
+++ b/src/cherab/imas/datasets/_builtin.py
@@ -91,7 +91,7 @@ def _iter_jintrac_radiation_values_data() -> tuple[IDSToplevel, IDSToplevel]:
     if subset_name is None:
         raise RuntimeError("Unable to find GGD subset id=5 (cells) in source grid_ggd.")
 
-    num_cells = len(subsets[subset_name])
+    num_cells = len(subsets[subset_name][0])
     if num_cells == 0:
         raise RuntimeError("The selected GGD subset (id=5) contains no cells.")
 

From 1f27c80472d2cbbc5a25f2fcd7865513c499c2a9 Mon Sep 17 00:00:00 2001
From: munechika-koyo <munechika.koyo@gmail.com>
Date: Tue, 11 Aug 2026 09:03:03 +0200
Subject: [PATCH 15/36] =?UTF-8?q?=F0=9F=8E=A8=20Update=20=5Fcreate=5Frad?=
 =?UTF-8?q?=5Ffunc=5Fggd:=20Enhance=20grid=20subset=20handling=20with=20va?=
 =?UTF-8?q?lidity=20mask?=
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit

---
 src/cherab/imas/emitter/radiation.py | 5 +++--
 1 file changed, 3 insertions(+), 2 deletions(-)

diff --git a/src/cherab/imas/emitter/radiation.py b/src/cherab/imas/emitter/radiation.py
index 2177818..1ffb3b4 100644
--- a/src/cherab/imas/emitter/radiation.py
+++ b/src/cherab/imas/emitter/radiation.py
@@ -144,8 +144,9 @@ def _create_rad_func_ggd(
     grid, subsets, subset_id = load_grid(grid_ggd, with_subsets=True)
     grid_subset_name, grid_subset_index = get_subset_name_index(subset_id, grid_subset_id)
 
-    if not np.array_equal(subsets[grid_subset_name], np.arange(grid.num_cell, dtype=int)):
-        grid = grid.subset(subsets[grid_subset_name], name=grid_subset_name)
+    subset_indices, subset_mask = subsets[grid_subset_name]
+    if not np.array_equal(subset_indices, np.arange(grid.num_cell, dtype=int)):
+        grid = grid.subset(subset_indices, name=grid_subset_name, valid_data_mask=subset_mask)
 
     rad_func = AxisymmetricMapper(grid.interpolator(data, **interp_kwargs))
 

From 588e2d036c52a17add9674789d4a2775e533f1f9 Mon Sep 17 00:00:00 2001
From: munechika-koyo <munechika.koyo@gmail.com>
Date: Tue, 11 Aug 2026 09:03:17 +0200
Subject: [PATCH 16/36] =?UTF-8?q?=E2=9C=A8=20Enhance=20species=20module:?=
 =?UTF-8?q?=20Add=20element=20styling,=20improve=20species=20type=20descri?=
 =?UTF-8?q?ptions,=20and=20implement=20new=20profile=20data=20methods?=
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit

---
 src/cherab/imas/ids/common/species.py | 310 +++++++++++++++++++++++---
 1 file changed, 280 insertions(+), 30 deletions(-)

diff --git a/src/cherab/imas/ids/common/species.py b/src/cherab/imas/ids/common/species.py
index 097c584..a9f2ca2 100644
--- a/src/cherab/imas/ids/common/species.py
+++ b/src/cherab/imas/ids/common/species.py
@@ -19,11 +19,16 @@
 
 from __future__ import annotations
 
-from dataclasses import dataclass, field
+from colorsys import hsv_to_rgb
+from dataclasses import dataclass, field, fields, is_dataclass
 from enum import Enum
+from typing import get_args, get_origin, get_type_hints
+from zlib import crc32
 
 import numpy as np
 from numpy.typing import NDArray
+from rich.text import Text
+from rich.tree import Tree
 
 from cherab.core.atomic.elements import Element, Isotope, lookup_isotope
 from imas.ids_defs import EMPTY_FLOAT, EMPTY_INT
@@ -36,6 +41,7 @@
     "ProfileData",
     "SpeciesComposition",
     "VelocityData",
+    "select_profile_data",
     "get_ion_state",
     "get_neutral_state",
     "get_ion",
@@ -44,6 +50,22 @@
 ]
 
 
+def _element_style(symbol: str) -> str:
+    """Return a stable bright style derived directly from an element symbol.
+
+    The hue excludes the red sector and does not depend on the size or ordering of a palette.
+
+    Returns
+    -------
+    str
+        Rich style containing a deterministic true-color foreground.
+    """
+    hash_fraction = (crc32(symbol.encode()) >> 20) / 0xFFF
+    hue = (35.0 + 290.0 * hash_fraction) / 360.0
+    rgb = tuple(round(channel * 255) for channel in hsv_to_rgb(hue, 0.6, 1.0))
+    return f"bold #{rgb[0]:02x}{rgb[1]:02x}{rgb[2]:02x}"
+
+
 class SpeciesType(Enum):
     """Enumeration of species types in IMAS."""
 
@@ -56,9 +78,9 @@ class SpeciesType(Enum):
     NEUTRAL_BUNDLE = "neutral_bundle"
     """Bundle of neutral states"""
     MOLECULE = "molecule"
-    """Single molecule state"""
+    """Single molecular state; neutral or charged according to ``z_min`` and ``z_max``"""
     MOLECULAR_BUNDLE = "molecular_bundle"
-    """Bundle of molecular states"""
+    """Bundle of molecular states or charge states"""
 
 
 @dataclass
@@ -69,9 +91,9 @@ class SpeciesData:
     """Minimum ionization state of the species"""
     z_max: int
     """Maximum ionization state of the species"""
-    element: Element | None = None
+    element: Element | Isotope | None = None
     """Element that makes up the species, if it is a single particle"""
-    elements: tuple[Element, ...] = field(default_factory=tuple)
+    elements: tuple[Element | Isotope, ...] = field(default_factory=tuple)
     """Elements that make up the species, if it is a molecule"""
     species_type: SpeciesType | None = None
     """Type of species"""
@@ -103,12 +125,93 @@ def __str__(self) -> str:
             else:
                 return f"{self.species_type.value} (z={self.z_min}-{self.z_max})"
         elif self.species_type == SpeciesType.MOLECULE:
-            return f"{'-'.join(el.symbol for el in self.elements)} {self.species_type.value}"
+            molecule = f"{'-'.join(el.symbol for el in self.elements)} {self.species_type.value}"
+            return f"{molecule} (z=+{self.z_min})" if self.z_min else molecule
         elif self.species_type == SpeciesType.MOLECULAR_BUNDLE:
             return f"{'-'.join(el.symbol for el in self.elements)} {self.species_type.value} (z={self.z_min}-{self.z_max})"
         else:
             return "Unknown species type"
 
+    def compact_label(self) -> str:
+        """Return a compact label for display beneath a species-type group.
+
+        Returns
+        -------
+        str
+            Elemental or molecular symbol followed by charge information when non-zero.
+        """
+        symbol = self._symbol_label()
+        if symbol is None:
+            return "Unknown"
+
+        if self.z_min != self.z_max:
+            minimum = self._charge_label(self.z_min, omit_unit=False)
+            maximum = self._charge_label(self.z_max, omit_unit=False)
+            return f"{symbol} {minimum}–{maximum}"
+        if self.z_min:
+            return f"{symbol} {self._charge_label(self.z_min)}"
+        return symbol
+
+    def _symbol_label(self) -> str | None:
+        """Return the elemental or molecular symbol used at the start of display labels.
+
+        Returns
+        -------
+        str or None
+            Elemental or molecular symbol, if available.
+        """
+        if self.element is not None:
+            return self.element.symbol
+        if self.elements:
+            return "-".join(element.symbol for element in self.elements)
+        return None
+
+    @staticmethod
+    def _charge_label(charge: int, *, omit_unit: bool = True) -> str:
+        """Format a charge number as a linearized ionic charge.
+
+        Returns
+        -------
+        str
+            Charge magnitude followed by its sign, omitting a unit magnitude when requested.
+        """
+        magnitude = abs(charge)
+        sign = "+" if charge >= 0 else "-"
+        return sign if omit_unit and magnitude == 1 else f"{magnitude}{sign}"
+
+    def rich_label(self, compact: bool = False) -> Text:
+        """Return a Rich label colored consistently by element symbol.
+
+        Parameters
+        ----------
+        compact
+            Use the compact context-aware label instead of the standalone description.
+
+        Returns
+        -------
+        `rich.text.Text`
+            Species label with an isotope-derived color when element data is available.
+        """
+        label = self.compact_label() if compact else str(self)
+        element = self.element or (self.elements[0] if self.elements else None)
+        if element is None:
+            return Text(label)
+        return Text(label, style=_element_style(element.symbol))
+
+    def __rich__(self) -> Text:
+        """Return a Rich label colored consistently by element symbol.
+
+        Charge states of the same element or isotope share a color, while different isotopes use
+        distinct symbols and therefore different colors. Homonuclear molecules use the color of
+        their constituent; heteronuclear molecules use the first constituent.
+
+        Returns
+        -------
+        `rich.text.Text`
+            Species label with an element-derived color when element data is available.
+        """
+        return self.rich_label()
+
 
 @dataclass
 class VelocityData:
@@ -145,6 +248,30 @@ class ProfileData:
     velocity: VelocityData | None = None
     """Bulk velocity data of the species."""
 
+    def array_shapes(self) -> tuple[tuple[str, tuple[int, ...]], ...]:
+        """Return paths and shapes for every array stored in this profile.
+
+        The dataclass hierarchy is traversed dynamically, so newly added array fields are
+        included without changing this method.
+
+        Returns
+        -------
+        tuple[tuple[str, tuple[int, ...]], ...]
+            Field paths and corresponding array shapes.
+        """
+        shapes: list[tuple[str, tuple[int, ...]]] = []
+
+        def collect(value: object, path: str = "") -> None:
+            if isinstance(value, np.ndarray):
+                shapes.append((path, value.shape))
+            elif is_dataclass(value) and not isinstance(value, type):
+                for data_field in fields(value):
+                    child_path = f"{path}.{data_field.name}" if path else data_field.name
+                    collect(getattr(value, data_field.name), child_path)
+
+        collect(self)
+        return tuple(shapes)
+
 
 @dataclass
 class SpeciesComposition:
@@ -161,10 +288,137 @@ class SpeciesComposition:
     neutral_bundle: list[ProfileData] = field(default_factory=list)
     """Neutral bundle profiles."""
     molecule: list[ProfileData] = field(default_factory=list)
-    """Molecule profiles."""
+    """Neutral and charged molecule profiles."""
     molecular_bundle: list[ProfileData] = field(default_factory=list)
     """Molecular bundle profiles."""
 
+    def _profile_groups(self) -> tuple[tuple[str, tuple[ProfileData, ...], bool], ...]:
+        """Return profile groups discovered from the dataclass type annotations.
+
+        Returns
+        -------
+        tuple[tuple[str, tuple[ProfileData, ...], bool], ...]
+            Field names, stored profiles, and whether the field is a profile collection.
+        """
+        type_hints = get_type_hints(type(self))
+        groups = []
+        for data_field in fields(self):
+            annotation = type_hints[data_field.name]
+            value = getattr(self, data_field.name)
+            is_profile = isinstance(annotation, type) and issubclass(annotation, ProfileData)
+            list_args = get_args(annotation) if get_origin(annotation) is list else ()
+            is_profile_list = (
+                len(list_args) == 1
+                and isinstance(list_args[0], type)
+                and issubclass(list_args[0], ProfileData)
+            )
+            if is_profile:
+                profiles = (value,)
+                is_collection = False
+            elif is_profile_list:
+                profiles = tuple(value)
+                is_collection = True
+            else:
+                continue
+            groups.append((data_field.name, profiles, is_collection))
+
+        return tuple(groups)
+
+    def __str__(self) -> str:
+        """Return a concise summary containing species names and profile shapes.
+
+        Returns
+        -------
+        str
+            Multiline summary of species names and array shapes grouped by type.
+        """
+        lines = [self.__class__.__name__]
+        for group_name, profiles, is_collection in self._profile_groups():
+            if not profiles:
+                continue
+            count = f" ({len(profiles)})" if len(profiles) > 1 else ""
+            lines.append(f"  {group_name}{count}")
+            for profile in profiles:
+                if is_collection:
+                    lines.append(f"    - {profile.species.compact_label()}")
+                indent = "      " if is_collection else "    "
+                lines.extend(
+                    f"{indent}{path}: shape={shape}" for path, shape in profile.array_shapes()
+                )
+
+        return "\n".join(lines)
+
+    def __rich__(self) -> Tree:
+        """Return a Rich tree containing species names and profile shapes.
+
+        Returns
+        -------
+        `rich.tree.Tree`
+            Tree representation rendered by Rich-enabled output.
+        """
+        tree = Tree(Text(self.__class__.__name__, style="bold cyan"))
+        for group_name, profiles, is_collection in self._profile_groups():
+            if not profiles:
+                continue
+            count = f" ({len(profiles)})" if len(profiles) > 1 else ""
+            branch = tree.add(Text(f"{group_name}{count}", style="bold"))
+            for profile in profiles:
+                parent = (
+                    branch.add(profile.species.rich_label(compact=True))
+                    if is_collection
+                    else branch
+                )
+                for path, shape in profile.array_shapes():
+                    parent.add(Text(f"{path}: shape={shape}", style="dim"))
+
+        return tree
+
+
+def select_profile_data(
+    composition: SpeciesComposition,
+    indices: NDArray[np.intp],
+    source_size: int,
+) -> None:
+    """Select valid source entries from all one-dimensional profile arrays.
+
+    Parameters
+    ----------
+    composition
+        Species composition whose profile arrays are updated in place.
+    indices
+        Positions of valid source entries in the selected grid subset.
+    source_size
+        Number of entries in the source grid subset before invalid cells are removed.
+    """
+    profiles = [composition.electron]
+    for group_name in (
+        "ion",
+        "ion_bundle",
+        "neutral",
+        "neutral_bundle",
+        "molecule",
+        "molecular_bundle",
+    ):
+        profiles.extend(getattr(composition, group_name))
+
+    def select_profile(profile: ProfileData) -> None:
+        for profile_field in fields(profile):
+            value = getattr(profile, profile_field.name)
+            if isinstance(value, np.ndarray) and value.ndim == 1 and value.size == source_size:
+                setattr(profile, profile_field.name, value[indices])
+            elif isinstance(value, VelocityData):
+                for velocity_field in fields(value):
+                    velocity = getattr(value, velocity_field.name)
+                    if (
+                        isinstance(velocity, np.ndarray)
+                        and velocity.ndim == 1
+                        and velocity.size == source_size
+                    ):
+                        setattr(value, velocity_field.name, velocity[indices])
+
+    for profile in profiles:
+        select_profile(profile)
+
 
 def get_ion_state(
     state: IDSStructure,
@@ -234,15 +488,17 @@ def get_ion_state(
 
     if len(elements) > 1:  # molecular ions and bundles
         species_data.elements = elements
-        if z_min == z_max == 0:
-            species_data.species_type = SpeciesType.NEUTRAL_BUNDLE
-        elif z_min == z_max:
+        if z_min == z_max:
             species_data.species_type = SpeciesType.MOLECULE
             species_data.vibrational_mode = (
-                str(state.vibrational_mode) if len(state.vibrational_mode) else None
+                str(getattr(state, "vibrational_mode", "")).strip()
+                if len(getattr(state, "vibrational_mode", ""))
+                else None
             )
             species_data.vibrational_level = (
-                state.vibrational_level if state.vibrational_level != EMPTY_FLOAT else None
+                getattr(state, "vibrational_level", EMPTY_FLOAT)
+                if getattr(state, "vibrational_level", EMPTY_FLOAT) != EMPTY_FLOAT
+                else None
             )
         else:
             species_data.species_type = SpeciesType.MOLECULAR_BUNDLE
@@ -281,25 +537,19 @@ def get_neutral_state(state: IDSStructure, elements: tuple[Element, ...]) -> Spe
         if len(getattr(state, "electron_configuration", "")) > 0
         else None,
     )
-    if len(elements) > 1:  # molecules and bundles
+    if len(elements) > 1:  # molecules
         species_data.elements = elements
-        if (
-            getattr(state, "vibrational_mode", None)
-            and getattr(state, "vibrational_level", EMPTY_FLOAT) != EMPTY_FLOAT
-        ):
-            species_data.species_type = SpeciesType.MOLECULE
-            species_data.vibrational_mode = (
-                str(getattr(state, "vibrational_mode", "")).strip()
-                if len(getattr(state, "vibrational_mode", ""))
-                else None
-            )
-            species_data.vibrational_level = (
-                getattr(state, "vibrational_level", EMPTY_FLOAT)
-                if getattr(state, "vibrational_level", EMPTY_FLOAT) != EMPTY_FLOAT
-                else None
-            )
-        else:
-            species_data.species_type = SpeciesType.NEUTRAL_BUNDLE
+        species_data.species_type = SpeciesType.MOLECULE
+        species_data.vibrational_mode = (
+            str(getattr(state, "vibrational_mode", "")).strip()
+            if len(getattr(state, "vibrational_mode", ""))
+            else None
+        )
+        species_data.vibrational_level = (
+            getattr(state, "vibrational_level", EMPTY_FLOAT)
+            if getattr(state, "vibrational_level", EMPTY_FLOAT) != EMPTY_FLOAT
+            else None
+        )
     else:  # neutrals
         species_data.element = elements[0]
         species_data.species_type = SpeciesType.NEUTRAL

From 7608570d31d7aa7c13b1574ec732d609b07042be Mon Sep 17 00:00:00 2001
From: munechika-koyo <munechika.koyo@gmail.com>
Date: Tue, 11 Aug 2026 09:04:38 +0200
Subject: [PATCH 17/36] =?UTF-8?q?=E2=9C=A8=20Enhance=20load=5Fprofiles:=20?=
 =?UTF-8?q?Add=20support=20for=20MOLECULE=20and=20MOLECULAR=5FBUNDLE=20spe?=
 =?UTF-8?q?cies=20types=20in=20core=20and=20edge=20species=20loading=20fun?=
 =?UTF-8?q?ctions?=
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit

---
 .../imas/ids/core_profiles/load_profiles.py   | 19 ++++++++++++++++---
 .../imas/ids/edge_profiles/load_profiles.py   | 19 ++++++++++++++++---
 2 files changed, 32 insertions(+), 6 deletions(-)

diff --git a/src/cherab/imas/ids/core_profiles/load_profiles.py b/src/cherab/imas/ids/core_profiles/load_profiles.py
index aef2dff..c92c0d9 100644
--- a/src/cherab/imas/ids/core_profiles/load_profiles.py
+++ b/src/cherab/imas/ids/core_profiles/load_profiles.py
@@ -195,6 +195,16 @@ def load_core_species(
                     composition.ion.append(profile_data)
                     ion_uuids.add(species_tuple)
 
+                # === Case: MOLECULE ===
+                elif species_data.species_type == SpeciesType.MOLECULE:
+                    composition.molecule.append(profile_data)
+                    ion_uuids.add(species_tuple)
+
+                # === Case: MOLECULAR_BUNDLE ===
+                elif species_data.species_type == SpeciesType.MOLECULAR_BUNDLE:
+                    composition.molecular_bundle.append(profile_data)
+                    ion_bundle_uuids.add(species_tuple)
+
                 # === Case: ION_BUNDLE ===
                 elif species_data.species_type == SpeciesType.ION_BUNDLE:
                     # Split ion bundles into individual ion states
@@ -277,14 +287,17 @@ def load_core_species(
                 if species_data.species_type == SpeciesType.ION:
                     composition.ion.append(profile_data)
                     ion_uuids.add(species_tuple)
+                elif species_data.species_type == SpeciesType.MOLECULE:
+                    composition.molecule.append(profile_data)
+                    ion_uuids.add(species_tuple)
                 else:
                     print(
                         f"Warning! Skipping non-bundled ion with unexpected species {species_data}"
                     )
 
-    # ----------------------------
-    # === Neutrals (molecules) ===
-    # ----------------------------
+    # ----------------
+    # === Neutrals ===
+    # ----------------
     for neutral in profile_1d.neutral:
         elements = get_elements(neutral.element)
         if not len(elements):
diff --git a/src/cherab/imas/ids/edge_profiles/load_profiles.py b/src/cherab/imas/ids/edge_profiles/load_profiles.py
index 38c312c..46d45d0 100644
--- a/src/cherab/imas/ids/edge_profiles/load_profiles.py
+++ b/src/cherab/imas/ids/edge_profiles/load_profiles.py
@@ -220,6 +220,16 @@ def load_edge_species(
                     composition.ion.append(profile_data)
                     ion_uuids.add(species_tuple)
 
+                # === Case: MOLECULE ===
+                elif species_data.species_type == SpeciesType.MOLECULE:
+                    composition.molecule.append(profile_data)
+                    ion_uuids.add(species_tuple)
+
+                # === Case: MOLECULAR_BUNDLE ===
+                elif species_data.species_type == SpeciesType.MOLECULAR_BUNDLE:
+                    composition.molecular_bundle.append(profile_data)
+                    ion_bundle_uuids.add(species_tuple)
+
                 # === Case: ION_BUNDLE ===
                 elif species_data.species_type == SpeciesType.ION_BUNDLE:
                     # Split ion bundles into individual ion states
@@ -302,14 +312,17 @@ def load_edge_species(
                 if species_data.species_type == SpeciesType.ION:
                     composition.ion.append(profile_data)
                     ion_uuids.add(species_tuple)
+                elif species_data.species_type == SpeciesType.MOLECULE:
+                    composition.molecule.append(profile_data)
+                    ion_uuids.add(species_tuple)
                 else:
                     print(
                         f"Warning! Skipping non-bundled ion with unexpected species {species_data}"
                     )
 
-    # ----------------------------
-    # === Neutrals (molecules) ===
-    # ----------------------------
+    # ----------------
+    # === Neutrals ===
+    # ----------------
     for neutral in ggd_struct.neutral:
         elements = get_elements(neutral.element)
         if not len(elements):

From be78363a4c858f044cbc80ac8d0daaa6482e9f41 Mon Sep 17 00:00:00 2001
From: munechika-koyo <munechika.koyo@gmail.com>
Date: Tue, 11 Aug 2026 09:05:04 +0200
Subject: [PATCH 18/36] =?UTF-8?q?=E2=9C=A8=20Enhance=20plasma=20loading:?=
 =?UTF-8?q?=20Refactor=20grid=20subset=20handling=20and=20add=20profile=20?=
 =?UTF-8?q?data=20selection=20for=20edge=20and=20core=20plasma=20functions?=
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit

---
 src/cherab/imas/plasma/blend.py   | 11 +++++--
 src/cherab/imas/plasma/edge.py    | 49 +++++++++++++++++++++++++++++--
 src/cherab/imas/plasma/utility.py | 20 ++++---------
 3 files changed, 59 insertions(+), 21 deletions(-)

diff --git a/src/cherab/imas/plasma/blend.py b/src/cherab/imas/plasma/blend.py
index 2489e2b..7a3abb0 100644
--- a/src/cherab/imas/plasma/blend.py
+++ b/src/cherab/imas/plasma/blend.py
@@ -37,11 +37,12 @@
 from ..ids.common import get_ids_time_slice
 from ..ids.common.ggd import load_grid
 from ..ids.common.grid_radial import get_psi_norm, load_core_grid
+from ..ids.common.species import select_profile_data
 from ..ids.core_profiles import load_core_species
 from ..ids.edge_profiles import load_edge_species
 from ..math.blend import blend_core_edge_functions
 from .core import get_core_interpolators, load_core_plasma
-from .edge import get_edge_interpolators, load_edge_plasma
+from .edge import _get_profile_indices, get_edge_interpolators, load_edge_plasma
 from .equilibrium import load_equilibrium, load_magnetic_field
 from .utility import (
     ZERO_VELOCITY,
@@ -284,9 +285,10 @@ def load_plasma(
     try:
         grid_subset_name, grid_subset_index = get_subset_name_index(subset_id, grid_subset_id)
 
-        if not np.array_equal(subsets[grid_subset_name], np.arange(grid.num_cell, dtype=int)):
+        subset_indices, subset_mask = subsets[grid_subset_name]
+        if not np.array_equal(subset_indices, np.arange(grid.num_cell, dtype=int)):
             # To reduce memory usage, create the sub-grid only if needed.
-            grid = grid.subset(subsets[grid_subset_name], name=grid_subset_name)
+            grid = grid.subset(subset_indices, name=grid_subset_name, valid_data_mask=subset_mask)
     except ValueError:
         print(
             f"Warning! Grid subset with identifier '{grid_subset_id}' not found in {subset_id}.",
@@ -299,6 +301,9 @@ def load_plasma(
         split_ion_bundles=split_ion_bundles,
         atomic_data=atomic_data,
     )
+    profile_indices, source_size = _get_profile_indices(grid_ggd, grid_subset_index)
+    if profile_indices is not None:
+        select_profile_data(composition_edge, profile_indices, source_size)
 
     # ----------------------------
     # === Create Plasma object ===
diff --git a/src/cherab/imas/plasma/edge.py b/src/cherab/imas/plasma/edge.py
index c4f8b7c..97d7603 100644
--- a/src/cherab/imas/plasma/edge.py
+++ b/src/cherab/imas/plasma/edge.py
@@ -39,7 +39,11 @@
 from ..ggd.base_mesh import GGDGrid
 from ..ids.common import get_ids_time_slice
 from ..ids.common.ggd import load_grid
-from ..ids.common.species import ProfileData, VelocityData
+from ..ids.common.species import (
+    ProfileData,
+    VelocityData,
+    select_profile_data,
+)
 from ..ids.edge_profiles import load_edge_species
 from ..math import UnitVector2D
 from .equilibrium import load_equilibrium, load_magnetic_field
@@ -160,9 +164,10 @@ def load_edge_plasma(
 
     try:
         grid_subset_name, grid_subset_index = get_subset_name_index(subset_id, grid_subset_id)
-        if not np.array_equal(subsets[grid_subset_name], np.arange(grid.num_cell, dtype=int)):
+        subset_indices, subset_mask = subsets[grid_subset_name]
+        if not np.array_equal(subset_indices, np.arange(grid.num_cell, dtype=int)):
             # To reduce memory usage, create the sub-grid only if needed.
-            grid = grid.subset(subsets[grid_subset_name], name=grid_subset_name)
+            grid = grid.subset(subset_indices, name=grid_subset_name, valid_data_mask=subset_mask)
     except ValueError:
         print(
             f"Warning! Grid subset with identifier '{grid_subset_id}' not found in {subset_id}.",
@@ -176,6 +181,9 @@ def load_edge_plasma(
         split_ion_bundles=split_ion_bundles,
         atomic_data=atomic_data,
     )
+    profile_indices, source_size = _get_profile_indices(grid_ggd_local, grid_subset_index)
+    if profile_indices is not None:
+        select_profile_data(composition, profile_indices, source_size)
 
     # ----------------------------
     # === Create Plasma object ===
@@ -276,6 +284,41 @@ def load_edge_plasma(
     return plasma
 
 
+def _get_profile_indices(
+    grid_ggd: IDSStructure, grid_subset_index: int
+) -> tuple[NDArray[np.intp] | None, int]:
+    """Return valid source positions for a selected GGD face subset.
+
+    Returns
+    -------
+    indices
+        Positions of valid face entries, or None when no GGD space is available.
+    source_size
+        Number of entries in the selected subset before filtering.
+    """
+    if not len(grid_ggd.space) or len(grid_ggd.space[0].objects_per_dimension) < 3:
+        return None, 0
+
+    faces = grid_ggd.space[0].objects_per_dimension[2].object
+    valid_faces = {
+        index
+        for index, face in enumerate(faces)
+        if face.has_value and face.nodes.has_value and len(face.nodes) >= 3
+    }
+    for subset in grid_ggd.grid_subset:
+        if subset.identifier.index.value != grid_subset_index:
+            continue
+        source_size = len(subset.element)
+        indices = [
+            position
+            for position, element in enumerate(subset.element)
+            if len(element.object) == 1 and element.object[0].index.value - 1 in valid_faces
+        ]
+        return np.asarray(indices, dtype=np.intp), source_size
+
+    return None, 0
+
+
 def get_edge_interpolators(
     grid: GGDGrid,
     profile: ProfileData,
diff --git a/src/cherab/imas/plasma/utility.py b/src/cherab/imas/plasma/utility.py
index 93f9f09..bb31e29 100644
--- a/src/cherab/imas/plasma/utility.py
+++ b/src/cherab/imas/plasma/utility.py
@@ -85,6 +85,10 @@ class ProfileInterpolator:
     """Interpolating function for the velocity profile."""
 
 
+# TODO: Add molecular species to Cherab Plasma when molecular models are available:
+# - represent neutral and charged molecules in Plasma.composition;
+# - construct distributions using molecular masses;
+# - integrate molecular collision and radiative models.
 def warn_unsupported_species(
     composition: SpeciesComposition,
     species_type: Literal["ion_bundle", "molecule", "molecular_bundle"],
@@ -100,21 +104,7 @@ def warn_unsupported_species(
     """
     profiles = getattr(composition, species_type, None)
     if profiles is not None and len(profiles) > 0:
-        names: list[str] = []
-        for profile_data in profiles:
-            name: str | None = getattr(profile_data.species, "name", None)
-            if name is None:
-                element = getattr(profile_data.species, "element", None)
-                if element is None:
-                    elements = getattr(profile_data.species, "elements", None)
-                    if elements is None or len(elements) == 0:
-                        name = "Unknown"
-                    else:
-                        name = "".join([e.name for e in elements])
-                else:
-                    name = element.name
-
-            names.append(name)
+        names = [str(profile_data.species) for profile_data in profiles]
 
         print(
             f"Warning! Species of type '{species_type}' are currently not supported.\n"

From 0e3dc72013d6d49ca955fc1f60f31b28403046d8 Mon Sep 17 00:00:00 2001
From: munechika-koyo <munechika.koyo@gmail.com>
Date: Tue, 11 Aug 2026 09:23:20 +0200
Subject: [PATCH 19/36] =?UTF-8?q?=E2=9C=A8=20Add=20comprehensive=20tests?=
 =?UTF-8?q?=20for=20species=20and=20core/edge=20profiles:=20Implement=20va?=
 =?UTF-8?q?lidation=20for=20molecular=20ions,=20cell=20data,=20and=20geome?=
 =?UTF-8?q?try=20calculations=20in=20new=20test=20files.?=
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit

---
 tests/ggd/test_interpolator_cache.py          |  12 +
 tests/ggd/test_load_unstruct_2d.py            |  36 +++
 tests/ggd/test_unstruct_2d_mesh.py            |  69 ++++++
 tests/ids/common/test_species.py              | 228 ++++++++++++++++++
 .../core_profiles/test_load_core_profiles.py  |  42 ++++
 .../edge_profiles/test_load_edge_profiles.py  |  21 ++
 tests/plasma/test_edge.py                     |  34 +++
 7 files changed, 442 insertions(+)
 create mode 100644 tests/ggd/test_load_unstruct_2d.py
 create mode 100644 tests/ggd/test_unstruct_2d_mesh.py
 create mode 100644 tests/ids/common/test_species.py
 create mode 100644 tests/ids/core_profiles/test_load_core_profiles.py
 create mode 100644 tests/ids/edge_profiles/test_load_edge_profiles.py

diff --git a/tests/ggd/test_interpolator_cache.py b/tests/ggd/test_interpolator_cache.py
index 2290302..a3c1ab1 100644
--- a/tests/ggd/test_interpolator_cache.py
+++ b/tests/ggd/test_interpolator_cache.py
@@ -12,6 +12,18 @@
 from cherab.imas.ggd.unstruct_3d_mesh import UnstructGrid3D
 
 
+def test_geometry_hash_supports_ragged_cells():
+    vertices = np.array([[0.0, 0.0], [1.0, 0.0], [1.0, 1.0], [0.0, 1.0], [2.0, 0.0]])
+    cells = [np.array([0, 1, 2, 3]), np.array([1, 4, 2])]
+
+    grid = UnstructGrid2D(vertices, cells, coordinate_system="cartesian")
+    same_grid = UnstructGrid2D(vertices.copy(), [cell.copy() for cell in cells])
+    changed_grid = UnstructGrid2D(vertices, [cells[0], np.array([1, 2, 4])])
+
+    assert grid._interpolator_geometry_hash() == same_grid._interpolator_geometry_hash()
+    assert grid._interpolator_geometry_hash() != changed_grid._interpolator_geometry_hash()
+
+
 def _assert_disk_cache_file_exists(
     cache_dir: Path,
     grid: UnstructGrid2D | UnstructGrid2DExtended | UnstructGrid3D,
diff --git a/tests/ggd/test_load_unstruct_2d.py b/tests/ggd/test_load_unstruct_2d.py
new file mode 100644
index 0000000..f67a6e7
--- /dev/null
+++ b/tests/ggd/test_load_unstruct_2d.py
@@ -0,0 +1,36 @@
+from typing import cast
+
+import numpy as np
+import pytest
+from imas import DBEntry
+
+from cherab.imas.ggd import CellData
+from cherab.imas.ggd.base_mesh import as_cell_data
+from cherab.imas.ids.common import get_ids_time_slice
+from cherab.imas.ids.common.ggd import load_grid
+
+
+def test_cell_data_validation():
+    """Cell data accepts one-dimensional float sequences and arrays only."""
+    data: CellData = (1.0, 2.0)
+
+    np.testing.assert_array_equal(as_cell_data(data, 2), np.array([1.0, 2.0]))
+    with pytest.raises(ValueError, match="one-dimensional"):
+        as_cell_data(cast(CellData, [[1.0, 2.0]]), 2)
+    with pytest.raises(ValueError, match="contain 2 values"):
+        as_cell_data(np.array([1.0]), 2)
+
+
+def test_iter_solps_cells_subset(path_iter_solps: str):
+    """The SOLPS ``Cells`` subset maps directly onto the loaded grid cells."""
+    with DBEntry(path_iter_solps, "r") as entry:
+        ids = get_ids_time_slice(entry, "edge_profiles", time=0)
+
+    grid, subsets, _ = load_grid(ids.grid_ggd[0], with_subsets=True)
+
+    cells, valid_mask = subsets["Cells"]
+    subset = grid.subset(cells, name="Cells", valid_data_mask=valid_mask)
+
+    assert len(cells) == grid.num_cell
+    assert any(len(cell) == 3 for cell in grid.cells)
+    assert subset.num_cell == len(cells)
diff --git a/tests/ggd/test_unstruct_2d_mesh.py b/tests/ggd/test_unstruct_2d_mesh.py
new file mode 100644
index 0000000..36b64bd
--- /dev/null
+++ b/tests/ggd/test_unstruct_2d_mesh.py
@@ -0,0 +1,69 @@
+import numpy as np
+import pytest
+
+from cherab.imas.ggd.unstruct_2d_mesh import UnstructGrid2D
+from cherab.imas.math.polygon import calculate_2d_cell_geometry
+
+
+def test_cylindrical_cell_geometry_uses_area_centroid_and_full_torus_volume():
+    # An asymmetric trapezoid whose area centroid differs from its vertex mean.
+    vertices = np.array([[1.0, 0.0], [3.0, 0.0], [2.0, 2.0], [1.0, 2.0]])
+    grid = UnstructGrid2D(vertices, [[0, 1, 2, 3]])
+
+    np.testing.assert_allclose(grid.cell_area, [3.0])
+    np.testing.assert_allclose(grid.cell_centre, [[16.0 / 9.0, 8.0 / 9.0]])
+    np.testing.assert_allclose(grid.cell_volume, [32.0 * np.pi / 3.0])
+
+    assert not grid.cell_area.flags["WRITEABLE"]
+    assert not grid.cell_centre.flags["WRITEABLE"]
+    assert not grid.cell_volume.flags["WRITEABLE"]
+
+
+def test_cylindrical_subset_preserves_cell_volumes():
+    vertices = np.array(
+        [
+            [1.0, 0.0],
+            [2.0, 0.0],
+            [2.0, 1.0],
+            [1.0, 1.0],
+            [3.0, 0.0],
+            [3.0, 1.0],
+        ]
+    )
+    grid = UnstructGrid2D(vertices, [[0, 1, 2, 3], [1, 4, 5, 2]])
+
+    subset = grid.subset([1])
+
+    np.testing.assert_allclose(subset.cell_volume, grid.cell_volume[[1]])
+    assert not subset.cell_volume.flags["WRITEABLE"]
+
+
+def test_cartesian_grid_does_not_calculate_cell_volumes():
+    vertices = np.array([[0.0, 0.0], [2.0, 0.0], [0.0, 1.0]])
+
+    grid = UnstructGrid2D(vertices, [[0, 1, 2]], coordinate_system="cartesian")
+
+    np.testing.assert_allclose(grid.cell_area, [1.0])
+    np.testing.assert_allclose(grid.cell_centre, [[2.0 / 3.0, 1.0 / 3.0]])
+    assert not hasattr(grid, "_cell_volume")
+
+
+def test_zero_area_cell_is_rejected():
+    vertices = np.array([[1.0, 0.0], [2.0, 0.0], [3.0, 0.0]])
+
+    with pytest.raises(ValueError, match="positive area"):
+        UnstructGrid2D(vertices, [[0, 1, 2]])
+
+
+def test_large_geometry_calculation_openmp_path():
+    vertices = np.array([[0.0, 0.0], [2.0, 0.0], [0.0, 1.0]])
+    triangles = np.array([[0, 1, 2]], dtype=np.int32)
+    cell_to_triangle = np.empty((100_000, 2), dtype=np.int32)
+    cell_to_triangle[:, 0] = 0
+    cell_to_triangle[:, 1] = 1
+
+    centres, areas = calculate_2d_cell_geometry(vertices, triangles, cell_to_triangle)
+
+    np.testing.assert_allclose(centres[:, 0], 2.0 / 3.0)
+    np.testing.assert_allclose(centres[:, 1], 1.0 / 3.0)
+    np.testing.assert_allclose(areas, 1.0)
diff --git a/tests/ids/common/test_species.py b/tests/ids/common/test_species.py
new file mode 100644
index 0000000..c255a91
--- /dev/null
+++ b/tests/ids/common/test_species.py
@@ -0,0 +1,228 @@
+from dataclasses import dataclass, field
+from types import SimpleNamespace
+from typing import cast
+
+import numpy as np
+import pytest
+from imas.ids_defs import EMPTY_FLOAT
+from imas.ids_structure import IDSStructure
+from rich.text import Text
+from rich.tree import Tree
+
+from cherab.core.atomic.elements import deuterium, helium, hydrogen, neon, tritium
+from cherab.imas.ids.common.species import (
+    ProfileData,
+    SpeciesComposition,
+    SpeciesData,
+    SpeciesType,
+    get_ion_state,
+    get_neutral_state,
+)
+
+
+@dataclass
+class ExtendedProfileData(ProfileData):
+    pressure: np.ndarray | None = None
+
+
+@dataclass
+class ExtendedSpeciesComposition(SpeciesComposition):
+    diagnostic: list[ExtendedProfileData] = field(default_factory=list)
+
+
+@pytest.mark.parametrize(
+    ("vibrational_mode", "vibrational_level"),
+    [
+        ("", EMPTY_FLOAT),
+        ("A_g", 1.0),
+    ],
+)
+def test_get_neutral_state_classifies_d2_as_molecule(
+    vibrational_mode: str,
+    vibrational_level: float,
+) -> None:
+    state = SimpleNamespace(
+        electron_configuration="",
+        vibrational_mode=vibrational_mode,
+        vibrational_level=vibrational_level,
+    )
+
+    species = get_neutral_state(cast(IDSStructure, state), (deuterium, deuterium))
+
+    assert species.species_type is SpeciesType.MOLECULE
+    assert species.elements == (deuterium, deuterium)
+    assert species.vibrational_mode == (vibrational_mode or None)
+    assert species.vibrational_level == (
+        vibrational_level if vibrational_level != EMPTY_FLOAT else None
+    )
+
+
+@pytest.mark.parametrize(
+    ("z_min", "z_max", "species_type"),
+    [
+        (0.0, 0.0, SpeciesType.MOLECULE),
+        (1.0, 1.0, SpeciesType.MOLECULE),
+        (1.0, 2.0, SpeciesType.MOLECULAR_BUNDLE),
+    ],
+)
+def test_get_ion_state_classifies_molecular_species(
+    z_min: float,
+    z_max: float,
+    species_type: SpeciesType,
+) -> None:
+    state = SimpleNamespace(
+        z_min=z_min,
+        z_max=z_max,
+        electron_configuration="",
+        vibrational_mode="",
+        vibrational_level=EMPTY_FLOAT,
+    )
+
+    species = get_ion_state(cast(IDSStructure, state), 0, (deuterium, deuterium))
+
+    assert species.species_type is species_type
+    assert species.elements == (deuterium, deuterium)
+    assert species.z_min == int(z_min)
+    assert species.z_max == int(z_max)
+
+
+def test_molecular_ion_string_includes_charge() -> None:
+    state = SimpleNamespace(
+        z_min=1.0,
+        z_max=1.0,
+        electron_configuration="",
+        vibrational_mode="",
+        vibrational_level=EMPTY_FLOAT,
+    )
+
+    species = get_ion_state(cast(IDSStructure, state), 0, (deuterium, deuterium))
+
+    assert str(species) == "D-D molecule (z=+1)"
+    assert species.compact_label() == "D-D +"
+
+
+@pytest.mark.parametrize(
+    ("z_min", "z_max", "expected"),
+    [
+        (0, 0, "Ne"),
+        (1, 1, "Ne +"),
+        (3, 3, "Ne 3+"),
+        (1, 10, "Ne 1+–10+"),
+    ],
+)
+def test_compact_label_uses_linearized_ionic_charge(z_min: int, z_max: int, expected: str) -> None:
+    species_type = SpeciesType.ION if z_min == z_max else SpeciesType.ION_BUNDLE
+    species = SpeciesData(
+        z_min=z_min,
+        z_max=z_max,
+        element=neon,
+        species_type=species_type,
+    )
+
+    assert species.compact_label() == expected
+
+
+def test_species_rich_color_is_shared_across_charge_and_molecular_states() -> None:
+    d_ion = SpeciesData(
+        z_min=1,
+        z_max=1,
+        element=deuterium,
+        species_type=SpeciesType.ION,
+    )
+    d2_ion = SpeciesData(
+        z_min=1,
+        z_max=1,
+        elements=(deuterium, deuterium),
+        species_type=SpeciesType.MOLECULE,
+    )
+    he_ion = SpeciesData(z_min=1, z_max=1, element=helium, species_type=SpeciesType.ION)
+
+    assert d_ion.__rich__().style == d2_ion.__rich__().style
+    assert d_ion.__rich__().style != he_ion.__rich__().style
+
+
+def test_species_rich_color_distinguishes_isotopes() -> None:
+    styles = {
+        SpeciesData(z_min=0, z_max=0, element=isotope, species_type=SpeciesType.NEUTRAL)
+        .__rich__()
+        .style
+        for isotope in (hydrogen, deuterium, tritium)
+    }
+
+    assert len(styles) == 3
+    assert all(str(style).startswith("bold #") for style in styles)
+
+
+def test_species_rich_color_includes_charge() -> None:
+    species = SpeciesData(z_min=3, z_max=3, element=neon, species_type=SpeciesType.ION)
+
+    label = species.rich_label(compact=True)
+
+    assert label.plain == "Ne 3+"
+    assert str(label.style).startswith("bold #")
+    assert not label.spans
+
+
+def test_species_composition_string_summarizes_species_and_profile_shapes() -> None:
+    density = np.ones(3)
+    electron = ProfileData(SpeciesData(z_min=-1, z_max=-1), density=density)
+    d2 = ProfileData(
+        SpeciesData(
+            z_min=0,
+            z_max=0,
+            elements=(deuterium, deuterium),
+            species_type=SpeciesType.MOLECULE,
+        ),
+        density=density,
+    )
+    d2_ion = ProfileData(
+        SpeciesData(
+            z_min=1,
+            z_max=1,
+            elements=(deuterium, deuterium),
+            species_type=SpeciesType.MOLECULE,
+        ),
+        density=density,
+    )
+    composition = SpeciesComposition(electron=electron, molecule=[d2, d2_ion])
+
+    expected = """SpeciesComposition
+  electron
+    density: shape=(3,)
+  molecule (2)
+    - D-D
+      density: shape=(3,)
+    - D-D +
+      density: shape=(3,)"""
+
+    assert str(composition) == expected
+    assert "[1. 1. 1.]" not in str(composition)
+    assert repr(composition).startswith("SpeciesComposition(electron=ProfileData(")
+
+    rich_tree = composition.__rich__()
+    assert isinstance(rich_tree, Tree)
+    assert isinstance(rich_tree.label, Text)
+    assert rich_tree.label.plain == "SpeciesComposition"
+    assert all(isinstance(branch.label, Text) for branch in rich_tree.children)
+    assert [cast(Text, branch.label).plain for branch in rich_tree.children] == [
+        "electron",
+        "molecule (2)",
+    ]
+    assert all(isinstance(node.label, Text) for node in rich_tree.children[1].children)
+    molecule_labels = [cast(Text, node.label).plain for node in rich_tree.children[1].children]
+    assert molecule_labels == ["D-D", "D-D +"]
+
+
+def test_species_summary_discovers_new_dataclass_fields() -> None:
+    electron = ProfileData(SpeciesData(z_min=-1, z_max=-1))
+    diagnostic = ExtendedProfileData(
+        SpeciesData(z_min=1, z_max=1, element=deuterium, species_type=SpeciesType.ION),
+        pressure=np.ones((2, 4)),
+    )
+    composition = ExtendedSpeciesComposition(electron=electron, diagnostic=[diagnostic])
+
+    summary = str(composition)
+
+    assert "  diagnostic\n" in summary
+    assert "diagnostic (1)" not in summary
+    assert "pressure: shape=(2, 4)" in summary
diff --git a/tests/ids/core_profiles/test_load_core_profiles.py b/tests/ids/core_profiles/test_load_core_profiles.py
new file mode 100644
index 0000000..d9a3691
--- /dev/null
+++ b/tests/ids/core_profiles/test_load_core_profiles.py
@@ -0,0 +1,42 @@
+from types import SimpleNamespace
+from typing import cast
+
+import numpy as np
+from imas.ids_defs import EMPTY_FLOAT
+from imas.ids_structure import IDSStructure
+
+from cherab.imas.ids.common.species import ProfileData, SpeciesType
+from cherab.imas.ids.core_profiles import load_profiles
+
+
+def test_load_core_species_preserves_molecular_ion(monkeypatch) -> None:
+    state = SimpleNamespace(
+        z_min=1.0,
+        z_max=1.0,
+        electron_configuration="",
+        vibrational_mode="",
+        vibrational_level=EMPTY_FLOAT,
+    )
+    molecular_ion = SimpleNamespace(
+        element=[SimpleNamespace(a=2.0, z_n=1.0, atoms_n=2)],
+        state=[state],
+    )
+    profile_1d = SimpleNamespace(electrons=object(), ion=[molecular_ion], neutral=[])
+
+    def fake_load_core_profiles(_structure, species, backup_species_struct=None):
+        return ProfileData(
+            species=species,
+            density=np.ones(1),
+            temperature=np.ones(1),
+        )
+
+    monkeypatch.setattr(load_profiles, "load_core_profiles", fake_load_core_profiles)
+    monkeypatch.setattr(load_profiles, "_get_profile", lambda *_args, **_kwargs: None)
+
+    composition = load_profiles.load_core_species(
+        cast(IDSStructure, profile_1d), split_ion_bundles=False
+    )
+
+    assert len(composition.molecule) == 1
+    assert composition.molecule[0].species.species_type is SpeciesType.MOLECULE
+    assert composition.molecule[0].species.z_min == 1
diff --git a/tests/ids/edge_profiles/test_load_edge_profiles.py b/tests/ids/edge_profiles/test_load_edge_profiles.py
new file mode 100644
index 0000000..ea6a025
--- /dev/null
+++ b/tests/ids/edge_profiles/test_load_edge_profiles.py
@@ -0,0 +1,21 @@
+from imas import DBEntry
+
+from cherab.imas.ids.common import get_ids_time_slice
+from cherab.imas.ids.common.species import SpeciesType
+from cherab.imas.ids.edge_profiles import load_edge_species
+
+
+def test_load_edge_species_preserves_solps_molecules(path_iter_solps: str) -> None:
+    with DBEntry(path_iter_solps, "r") as entry:
+        ids = get_ids_time_slice(entry, "edge_profiles", time=0)
+
+    composition = load_edge_species(ids.ggd[0], split_ion_bundles=False)
+    d2_charge_states = {
+        profile.species.z_min
+        for profile in composition.molecule
+        if profile.species.species_type is SpeciesType.MOLECULE
+        and len(profile.species.elements) == 2
+    }
+
+    assert d2_charge_states == {0, 1}
+    assert not composition.neutral_bundle
diff --git a/tests/plasma/test_edge.py b/tests/plasma/test_edge.py
index 964ee1e..51998be 100644
--- a/tests/plasma/test_edge.py
+++ b/tests/plasma/test_edge.py
@@ -49,6 +49,40 @@ def test_load_edge_plasma(path_iter_jintrac: str):
     assert len(ion_charges) == neon.atomic_number + 1
 
 
+@pytest.mark.parametrize(
+    "grid_subset_id",
+    [
+        pytest.param(5, id="Cells"),
+        pytest.param(-1, id="Inner-core"),
+        pytest.param(-2, id="Inner-SOL"),
+        pytest.param(-3, id="Lower-inner-divertor"),
+        pytest.param(-5, id="Outer-core"),
+        pytest.param(-6, id="Outer-SOL"),
+        pytest.param(-8, id="Lower-outer-divertor"),
+        pytest.param(-101, id="Neutral-pressure-cells"),
+        pytest.param(22, id="CORE"),
+        pytest.param(23, id="SOL"),
+        pytest.param(24, id="OUTER_DIVERTOR"),
+        pytest.param(25, id="INNER_DIVERTOR"),
+        pytest.param(12, id="Inner-Midplane"),
+        pytest.param(11, id="Outer-Midplane"),
+    ],
+)
+def test_load_edge_plasma_solps(path_iter_solps: str, grid_subset_id: int):
+    """Test loading edge plasma from the SOLPS dataset."""
+    plasma = load_edge_plasma(
+        path_iter_solps,
+        "r",
+        grid_subset_id=grid_subset_id,
+        split_ion_bundles=False,
+    )
+
+    assert isinstance(plasma, Plasma)
+    assert plasma.geometry is not None
+    assert plasma.electron_distribution is not None
+    assert len(plasma.composition) > 0
+
+
 def test_load_edge_plasma_with_time(path_iter_jintrac: str):
     """Test loading edge plasma with specific time parameter."""
     plasma = load_edge_plasma(path_iter_jintrac, "r", time=0.1, split_ion_bundles=False)

From 4f88e5a81386cb2c1d39af9db7b10015693ab47d Mon Sep 17 00:00:00 2001
From: munechika-koyo <munechika.koyo@gmail.com>
Date: Tue, 11 Aug 2026 09:25:12 +0200
Subject: [PATCH 20/36] Edit edge plasma notebook to use ITER-SOLPS sample data

---
 docs/notebooks/plasma/2_edge_plasma.ipynb | 183 +++++++++++++++++-----
 1 file changed, 144 insertions(+), 39 deletions(-)

diff --git a/docs/notebooks/plasma/2_edge_plasma.ipynb b/docs/notebooks/plasma/2_edge_plasma.ipynb
index d79e8ba..d8c069f 100644
--- a/docs/notebooks/plasma/2_edge_plasma.ipynb
+++ b/docs/notebooks/plasma/2_edge_plasma.ipynb
@@ -10,7 +10,7 @@
     "This notebook demonstrates how to load and visualize edge plasma profiles using the `cherab.imas` interface.\n",
     "Here, we propose how to visualize edge plasmas with grid meshes defined in the IMAS data structure.\n",
     "\n",
-    "The example test data was calculated by JINTRAC for an ITER 15 MA H-mode scenario."
+    "The example test data was calculated by SOLPS-ITER for an ITER 15 MA H-mode scenario."
    ]
   },
   {
@@ -26,7 +26,7 @@
     "from matplotlib.colors import SymLogNorm\n",
     "from rich import print as rprint\n",
     "\n",
-    "from cherab.imas.datasets import iter_jintrac\n",
+    "from cherab.imas.datasets import iter_jintrac, iter_solps\n",
     "from cherab.imas.ggd import GGDGrid\n",
     "from cherab.imas.ids.common import get_ids_time_slice\n",
     "from cherab.imas.ids.common.ggd import load_grid\n",
@@ -115,7 +115,7 @@
    "id": "4",
    "metadata": {},
    "source": [
-    "## Retrieve ITER JINTRAC sample data"
+    "## Retrieve the sample data"
    ]
   },
   {
@@ -125,7 +125,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "path = iter_jintrac()"
+    "path = iter_solps()"
    ]
   },
   {
@@ -135,10 +135,10 @@
    "source": [
     "## Load grid and species data\n",
     "\n",
-    "### Select grid subset\n",
+    "### Plot all grid subsets\n",
     "\n",
     "In [edge_profiles IDS](https://imas-data-dictionary.readthedocs.io/en/latest/generated/ids/edge_profiles.html), there are multiple grid subsets defined.\n",
-    "Here, we choose the `\"cells\"` subset to visualize the edge plasma profiles."
+    "Here, we see what grid subsets are available and plot them all."
    ]
   },
   {
@@ -163,32 +163,64 @@
     ")\n",
     "\n",
     "# Print available grid subsets\n",
-    "rprint(\"Available grid subsets:\", subset_id)\n",
-    "\n",
-    "# Extract only \"cells\" subset\n",
-    "grid = grid.subset(subsets[\"cells\"])"
+    "rprint(\"Available grid subsets:\", subset_id)"
    ]
   },
   {
-   "cell_type": "markdown",
+   "cell_type": "code",
+   "execution_count": null,
    "id": "8",
    "metadata": {},
+   "outputs": [],
    "source": [
-    "### Load edge species data"
+    "subset_groups = [\n",
+    "    [\n",
+    "        \"Cells\",\n",
+    "    ],\n",
+    "    [\n",
+    "        \"Inner core\",\n",
+    "        \"Outer core\",\n",
+    "        \"Inner SOL\",\n",
+    "        \"Outer SOL\",\n",
+    "        \"Lower inner divertor\",\n",
+    "        \"Lower outer divertor\",\n",
+    "    ],\n",
+    "    [\n",
+    "        \"CORE\",\n",
+    "        \"SOL\",\n",
+    "        \"OUTER_DIVERTOR\",\n",
+    "        \"INNER_DIVERTOR\",\n",
+    "        \"Inner Midplane\",\n",
+    "        \"Outer Midplane\",\n",
+    "        \"Neutral pressure cells\",\n",
+    "    ],\n",
+    "]\n",
+    "\n",
+    "fig, axs = uplt.subplots(ncols=len(subset_groups))\n",
+    "\n",
+    "for ax, subset_names in zip(axs, subset_groups, strict=True):\n",
+    "    for i, subset_name in enumerate(subset_names):\n",
+    "        grid_subset = grid.subset(subsets[subset_name])\n",
+    "        grid_subset.plot_mesh(ax=ax, label=subset_name, edgecolor=f\"C{i}\")\n",
+    "\n",
+    "    ax.legend(ncols=1, loc=\"center\")\n",
+    "\n",
+    "axs.format(\n",
+    "    xlim=(4.0, 8.5),\n",
+    "    ylim=(-4.7, 4.8),\n",
+    "    grid=True,\n",
+    "    xlocator=1,\n",
+    "    ylocator=1,\n",
+    "    tickminor=True,\n",
+    ")"
    ]
   },
   {
-   "cell_type": "code",
-   "execution_count": null,
+   "cell_type": "markdown",
    "id": "9",
    "metadata": {},
-   "outputs": [],
    "source": [
-    "composition = load_edge_species(\n",
-    "    ids.ggd[0],\n",
-    "    grid_subset_index=subset_id[\"cells\"],\n",
-    "    split_ion_bundles=False,\n",
-    ")"
+    "### Load edge species data"
    ]
   },
   {
@@ -196,26 +228,33 @@
    "id": "10",
    "metadata": {},
    "source": [
-    "## Plot edge plasma profiles"
+    "We choose the `\"Cells\"` subset covering the entire edge region and load the corresponding edge species data. The `edge_profiles` IDS contains multiple species, and we can choose which one to visualize."
    ]
   },
   {
-   "cell_type": "markdown",
+   "cell_type": "code",
+   "execution_count": null,
    "id": "11",
    "metadata": {},
+   "outputs": [],
    "source": [
-    "### Grid mesh"
+    "grid_cells = grid.subset(subsets[\"Cells\"])\n",
+    "\n",
+    "composition = load_edge_species(\n",
+    "    ids.ggd[0],\n",
+    "    grid_subset_index=subset_id[\"Cells\"],\n",
+    "    split_ion_bundles=False,\n",
+    ")\n",
+    "\n",
+    "rprint(composition)"
    ]
   },
   {
-   "cell_type": "code",
-   "execution_count": null,
+   "cell_type": "markdown",
    "id": "12",
    "metadata": {},
-   "outputs": [],
    "source": [
-    "fig, ax = uplt.subplots()\n",
-    "ax = grid.plot_mesh(ax=ax)"
+    "## Plot edge plasma profiles"
    ]
   },
   {
@@ -241,20 +280,34 @@
     "fig, ax = uplt.subplots()\n",
     "ax = plot_grid_quantity(\n",
     "    ax,\n",
-    "    grid,\n",
+    "    grid_cells,\n",
     "    composition.electron.density,\n",
     "    title_center=\"Electron density\\n$n_\\\\mathrm{e}$ [m$^{-3}$]\",\n",
     "    logscale=True,\n",
     ")\n",
+    "ax.format(\n",
+    "    xlim=(4.0, 8.5),\n",
+    "    ylim=(-4.7, 4.8),\n",
+    "    xlocator=1,\n",
+    "    ylocator=1,\n",
+    "    tickminor=True,\n",
+    ")\n",
     "\n",
     "# Electron temperature\n",
     "fig, ax = uplt.subplots()\n",
     "ax = plot_grid_quantity(\n",
     "    ax,\n",
-    "    grid,\n",
+    "    grid_cells,\n",
     "    composition.electron.temperature,\n",
     "    title_center=\"Electron temperature\\n$T_\\\\mathrm{e}$ [eV]\",\n",
     "    logscale=True,\n",
+    ")\n",
+    "ax.format(\n",
+    "    xlim=(4.0, 8.5),\n",
+    "    ylim=(-4.7, 4.8),\n",
+    "    xlocator=1,\n",
+    "    ylocator=1,\n",
+    "    tickminor=True,\n",
     ")"
    ]
   },
@@ -278,13 +331,19 @@
     "\n",
     "for profile in composition.ion + composition.neutral + composition.molecule:\n",
     "    charge = profile.species.z_min\n",
-    "    element = profile.species.element or profile.species.elements[0]\n",
+    "    if (element := profile.species.element) is not None:\n",
+    "        symbol = element.symbol\n",
+    "    elif profile.species.elements:\n",
+    "        symbol = \"-\".join(element.symbol for element in profile.species.elements)\n",
+    "    else:\n",
+    "        symbol = \"Unknown\"\n",
+    "\n",
     "    if charge == 0:\n",
-    "        name = element.symbol\n",
+    "        name = symbol\n",
     "    elif charge == 1:\n",
-    "        name = f\"{element.symbol}$^+$\"\n",
+    "        name = f\"{symbol}$^+$\"\n",
     "    else:\n",
-    "        name = f\"{element.symbol}$^{{{charge}+}}$\"\n",
+    "        name = f\"{symbol}$^{{{charge}+}}$\"\n",
     "\n",
     "    # Density\n",
     "    data.append(\n",
@@ -296,7 +355,7 @@
     "            ),\n",
     "        )\n",
     "    )\n",
-    "    if element.atomic_number == 1:\n",
+    "    if (element := profile.species.element) is not None and element.atomic_number == 1:\n",
     "        # Temperature\n",
     "        if profile.temperature is not None and np.any(profile.temperature):\n",
     "            data.append(\n",
@@ -344,7 +403,7 @@
     "for i_ax, (quantity, kwargs) in enumerate(data):\n",
     "    ax = plot_grid_quantity(\n",
     "        axes[i_ax],\n",
-    "        grid,\n",
+    "        grid_cells,\n",
     "        quantity,\n",
     "        **kwargs,\n",
     "        cbar_kwargs=dict(\n",
@@ -362,6 +421,7 @@
     "    xtickloc=\"neither\",\n",
     "    ytickloc=\"neither\",\n",
     "    linestyle=\"none\",\n",
+    "    grid=False,\n",
     ")"
    ]
   },
@@ -384,10 +444,55 @@
     "Here, we demonstrate how to split the bundled species profiles into each charge state using the :obj:`.solve_coronal_equilibrium` function."
    ]
   },
+  {
+   "cell_type": "markdown",
+   "id": "19",
+   "metadata": {},
+   "source": [
+    "### Retrieve the bundled species profiles"
+   ]
+  },
   {
    "cell_type": "code",
    "execution_count": null,
-   "id": "19",
+   "id": "20",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "path = iter_jintrac()\n",
+    "\n",
+    "# Load edge_profiles IDs\n",
+    "with DBEntry(path, \"r\") as entry:\n",
+    "    ids = get_ids_time_slice(\n",
+    "        entry,\n",
+    "        \"edge_profiles\",\n",
+    "        time=0,\n",
+    "    )\n",
+    "\n",
+    "# Load grid object\n",
+    "grid, subsets, subset_id = load_grid(\n",
+    "    ids.grid_ggd[0],\n",
+    "    with_subsets=True,\n",
+    ")\n",
+    "\n",
+    "# Print available grid subsets\n",
+    "rprint(\"Available grid subsets:\", subset_id)\n",
+    "\n",
+    "grid_cells = grid.subset(subsets[\"cells\"])\n",
+    "\n",
+    "composition = load_edge_species(\n",
+    "    ids.ggd[0],\n",
+    "    grid_subset_index=subset_id[\"cells\"],\n",
+    "    split_ion_bundles=False,\n",
+    ")\n",
+    "\n",
+    "rprint(composition)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "21",
    "metadata": {},
    "outputs": [],
    "source": [
@@ -418,7 +523,7 @@
     "):\n",
     "    ax = plot_grid_quantity(\n",
     "        axes[i_ax],\n",
-    "        grid,\n",
+    "        grid_cells,\n",
     "        densities[i_ax, :],\n",
     "        title_center=f\"{neon.symbol}$^{{{charge}+}}$ density [m$^{{-3}}$]\",\n",
     "        logscale=True,\n",
@@ -457,7 +562,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.14.3"
+   "version": "3.14.6"
   }
  },
  "nbformat": 4,

From dfc8c40cfe5717a5ca7711d5f5b36bbe95472c53 Mon Sep 17 00:00:00 2001
From: munechika-koyo <munechika.koyo@gmail.com>
Date: Tue, 11 Aug 2026 09:41:00 +0200
Subject: [PATCH 21/36] =?UTF-8?q?=F0=9F=93=9D=20Update=20changelog?=
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit

---
 CHANGELOG.md | 29 ++++++++++++++++++++++++++++-
 1 file changed, 28 insertions(+), 1 deletion(-)

diff --git a/CHANGELOG.md b/CHANGELOG.md
index de1d9f8..51a012f 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -5,7 +5,7 @@ All notable changes to this project will be documented in this file.
 The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.1.0/),
 and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0.html).
 
-## [0.5.1]
+## [0.6.0]
 
 ### Added
 
@@ -13,18 +13,45 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 - Extend radiation emitter loading to support dual emissivity sources with improved validation
 - Add IDS path utility helpers and `get_entry_reference` for resolving entry references
 - Add unit tests for IDS path handling, grid loading via path references, and plasma utility entry-reference workflows
+- Add public `CellConnectivity`, `CellData`, and `VertexIndices` types and the `as_cell_data()`
+  validation helper for GGD meshes
+- Add validity-mask support for retaining the relationship between compacted 2D grids and their
+  source GGD face data
+- Add an optimized Cython implementation for calculating polygonal cell areas and area-weighted
+  centroids, with OpenMP support for large meshes
+- Add profile-shape inspection and plain-text/Rich summaries for species compositions
+- Add regression tests for 2D grid geometry, subset validity, molecular species, and core/edge
+  profile loading
+- Add a 2D radiation-emitter example notebook
 
 ### Changed
 
 - Improve radiation emitter loading checks for duplicate emissivity values and core-profile grid data
 - Refactor `load_grid` to support explicit `entry` selection and improve referenced-grid source resolution
 - Enhance plasma and emitter loading paths to use entry-reference aware grid resolution
+- **Breaking:** Change `load_grid(..., with_subsets=True)` and
+  `load_unstruct_grid_2d(..., with_subsets=True)` subset values from index arrays to
+  `(indices, valid_data_mask)` tuples
+- Extend `UnstructGrid2D` with source-data validity tracking so interpolation and plotting accept
+  either compacted cell data or source-sized data
+- Calculate cylindrical 2D cell volumes from area-weighted centroids over a full toroidal rotation
+- Filter edge and blended-plasma profile arrays consistently when invalid GGD faces are omitted
+- Classify neutral and charged molecular species consistently and retain molecular bundles in core
+  and edge species compositions
 - Promote `rich` from a test-only dependency to a runtime dependency
+- Summarize non-empty species groups with compact labels, conventional ionic charge notation, and
+  stable symbol-derived coloring in plain and `rich` tree output
+- Update the edge-plasma example to use the ITER-SOLPS sample dataset
 
 ### Fixed
 
 - Fix grid data loading checks for radiation core profiles
 - Improve error handling in radiation emitter loading workflows
+- Ignore missing, incomplete, and out-of-range GGD faces while preserving correct subset mappings
+- Fix triangle indexing, cylindrical cell-volume calculation, and volume preservation in 2D grid
+  subsets
+- Fix the Cython compile and link arguments on macOS
+- Preserve neutral and charged molecular species when loading core and edge compositions
 
 ## [0.5.0] - 2026-06-24
 

From 43bf231ea369b42a3100083c545ea1bb2babd9dd Mon Sep 17 00:00:00 2001
From: munechika-koyo <munechika.koyo@gmail.com>
Date: Tue, 11 Aug 2026 11:14:49 +0200
Subject: [PATCH 22/36] =?UTF-8?q?=E2=9C=A8=20Add=20internal=20helpers=20fo?=
 =?UTF-8?q?r=20opening=20IMAS=20database=20entries:=20Implement=20=5Fopen?=
 =?UTF-8?q?=5Fdbentry=5Ffor=5Freading=20and=20=5Fvalidate=5Fread=5Fmode=20?=
 =?UTF-8?q?functions=20for=20improved=20data=20loading.?=
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit

---
 src/cherab/imas/_dbentry.py | 68 +++++++++++++++++++++++++++++++++++++
 1 file changed, 68 insertions(+)
 create mode 100644 src/cherab/imas/_dbentry.py

diff --git a/src/cherab/imas/_dbentry.py b/src/cherab/imas/_dbentry.py
new file mode 100644
index 0000000..14ff621
--- /dev/null
+++ b/src/cherab/imas/_dbentry.py
@@ -0,0 +1,68 @@
+"""Internal helpers for opening IMAS database entries."""
+
+from __future__ import annotations
+
+from os import PathLike, fspath
+from typing import Any
+from warnings import warn
+
+from imas import DBEntry  # type: ignore[attr-defined]
+
+
+def _open_dbentry_for_reading(*args: Any, **kwargs: Any) -> DBEntry:
+    """Create a DBEntry for use by a CHERAB-IMAS data loader.
+
+    URI-style DBEntry construction requires an explicit mode in IMAS-Python. CHERAB-IMAS
+    loader APIs only read existing data, so this helper supplies ``"r"`` automatically.
+    The legacy ``(backend_id, db_name, pulse, run, ...)`` constructor is passed through
+    unchanged because it opens an existing entry when used as a context manager.
+
+    An explicitly supplied ``"r"`` is accepted temporarily for backwards compatibility.
+    Modes that may create or replace data are rejected.
+
+    Returns
+    -------
+    `imas.DBEntry`
+        Database entry opened or prepared for reading.
+
+    Raises
+    ------
+    TypeError
+        If mode is supplied both positionally and by keyword.
+    """
+    positional = list(args)
+
+    if positional and isinstance(positional[0], (str, PathLike)):
+        uri = fspath(positional.pop(0))
+        positional_mode = positional.pop(0) if positional else None
+        keyword_mode = kwargs.pop("mode", None)
+
+        if positional_mode is not None and keyword_mode is not None:
+            raise TypeError("DBEntry mode was specified both positionally and by keyword.")
+
+        mode = positional_mode if positional_mode is not None else keyword_mode
+        _validate_read_mode(mode)
+        return DBEntry(uri, "r", *positional, **kwargs)
+
+    if not positional and "uri" in kwargs:
+        uri = fspath(kwargs.pop("uri"))
+        mode = kwargs.pop("mode", None)
+        _validate_read_mode(mode)
+        return DBEntry(uri, "r", **kwargs)
+
+    return DBEntry(*positional, **kwargs)
+
+
+def _validate_read_mode(mode: Any) -> None:
+    if mode is None:
+        return
+
+    if mode != "r":
+        raise ValueError(f"CHERAB-IMAS loader APIs only support mode 'r'; received mode {mode!r}.")
+
+    warn(
+        "Passing mode 'r' to CHERAB-IMAS loader APIs is deprecated; "
+        "the read mode is now selected automatically.",
+        DeprecationWarning,
+        stacklevel=3,
+    )

From a8417b72508c3b75846ca1e8af683f23516b2577 Mon Sep 17 00:00:00 2001
From: munechika-koyo <munechika.koyo@gmail.com>
Date: Tue, 11 Aug 2026 11:29:30 +0200
Subject: [PATCH 23/36] =?UTF-8?q?=E2=99=BB=EF=B8=8F=20Refactor=20database?=
 =?UTF-8?q?=20entry=20handling:=20Replace=20DBEntry=20with=20`=5Fopen=5Fdb?=
 =?UTF-8?q?entry=5Ffor=5Freading`?=
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit

---
 src/cherab/imas/emitter/radiation.py  | 22 ++++++++++++----------
 src/cherab/imas/observer/bolometer.py | 18 ++++++++++--------
 src/cherab/imas/plasma/blend.py       | 24 ++++++++++++++----------
 src/cherab/imas/plasma/core.py        | 10 ++++++----
 src/cherab/imas/plasma/edge.py        | 10 ++++++----
 src/cherab/imas/plasma/equilibrium.py | 18 +++++++++++-------
 src/cherab/imas/wall/wall.py          | 23 +++++++++++++----------
 7 files changed, 72 insertions(+), 53 deletions(-)

diff --git a/src/cherab/imas/emitter/radiation.py b/src/cherab/imas/emitter/radiation.py
index 1ffb3b4..49208b4 100644
--- a/src/cherab/imas/emitter/radiation.py
+++ b/src/cherab/imas/emitter/radiation.py
@@ -35,10 +35,10 @@
 from cherab.core.math import AxisymmetricMapper
 from cherab.tools.emitters import RadiationFunction
 from cherab.tools.equilibrium import EFITEquilibrium
-from imas import DBEntry
 from imas.ids_struct_array import IDSStructArray
 from imas.ids_structure import IDSStructure
 
+from .._dbentry import _open_dbentry_for_reading
 from ..ggd import UnstructGrid2DExtended
 from ..ggd.base_mesh import InterpolatorCacheMode
 from ..ids.common import get_ids_time_slice, load_ids_path_reference, resolve_ids_path_reference
@@ -174,7 +174,7 @@ def _resolve_grid_ggd_reference(
             raise RuntimeError(
                 "Unable to resolve external grid_ggd.path without DBEntry arguments."
             )
-        with DBEntry(*db_args, **(db_kwargs or {})) as entry:
+        with _open_dbentry_for_reading(*db_args, **(db_kwargs or {})) as entry:
             resolved = load_ids_path_reference(entry, path)
     else:
         ids_name = ids_root.metadata.name
@@ -244,17 +244,19 @@ def load_radiation_emitter(
     Parameters
     ----------
     *args
-        Positional arguments passed to `imas.DBEntry`.
+        IMAS URI, netCDF path, or legacy positional arguments for `imas.DBEntry`.
+        For a URI or path, read mode is selected automatically; do not pass ``"r"``.
     time
         Time slice to load from the IDS, by default 0.0.
     occurrence
         Occurrence of the radiation IDS, by default 0.
     args2
-        Arguments passed to `imas.DBEntry` for the second emissivity. If None, the second emissivity
-        is not loaded, by default None.
+        URI, netCDF path, or legacy positional DBEntry arguments for the second emissivity.
+        Read mode is selected automatically. If None, the second emissivity is not loaded,
+        by default None.
     kwargs2
-        Keyword arguments passed to `imas.DBEntry` for the second emissivity. If None, the second
-        emissivity is not loaded, by default None.
+        Additional DBEntry options for the second emissivity. If None, no options are used,
+        by default None.
     time2
         Time slice to load for the second emissivity. By default, uses the same time as the first
         emissivity.
@@ -305,7 +307,7 @@ def load_radiation_emitter(
         Directory used when ``interpolator_cache="disk"``, by default None
         (uses the system cache directory, e.g., ``~/.cache/cherab/imas/interpolators``).
     **kwargs
-        Additional keyword arguments passed to `imas.DBEntry`.
+        Additional `imas.DBEntry` options, such as ``dd_version`` or ``xml_path``.
 
     Returns
     -------
@@ -346,7 +348,7 @@ def load_radiation_emitter(
     zmin: float = 0.0
 
     try:
-        with DBEntry(*args, **kwargs) as entry:
+        with _open_dbentry_for_reading(*args, **kwargs) as entry:
             ids = get_ids_time_slice(
                 entry,
                 "radiation",
@@ -360,7 +362,7 @@ def load_radiation_emitter(
 
     if args2 is not None and source != "coefficients":
         try:
-            with DBEntry(*args2, **(kwargs2 or {})) as entry:
+            with _open_dbentry_for_reading(*args2, **(kwargs2 or {})) as entry:
                 ids2 = get_ids_time_slice(
                     entry,
                     "radiation",
diff --git a/src/cherab/imas/observer/bolometer.py b/src/cherab/imas/observer/bolometer.py
index e68290a..aa71bd0 100644
--- a/src/cherab/imas/observer/bolometer.py
+++ b/src/cherab/imas/observer/bolometer.py
@@ -33,8 +33,8 @@
 from raysect.primitive.csg import CSGPrimitive
 
 from cherab.tools.observers.bolometry import BolometerCamera, BolometerFoil, BolometerSlit
-from imas.db_entry import DBEntry
 
+from .._dbentry import _open_dbentry_for_reading
 from ..ids.bolometer import load_cameras
 from ..ids.bolometer._camera import BoloCamera, Geometry
 from ..ids.bolometer.utility import CameraType, GeometryType
@@ -69,7 +69,7 @@ def load_bolometers(
 @overload
 def load_bolometers(
     uri: str,
-    mode: str,
+    mode: Literal["r"] | None = None,
     *,
     parent: _NodeBase | None = None,
     dd_version: str | None = None,
@@ -90,12 +90,14 @@ def load_bolometers(
     Parameters
     ----------
     *args
-        Arguments passed to `~imas.db_entry.DBEntry`.
+        IMAS URI, netCDF path, or legacy positional arguments for
+        `~imas.db_entry.DBEntry`. For a URI or path, read mode is selected automatically;
+        do not pass ``"r"``.
     parent
         The parent node of `~cherab.tools.observers.bolometry.BolometerCamera` in the Raysect
         scene-graph, by default None.
     **kwargs
-        Keyword arguments passed to `~imas.db_entry.DBEntry` constructor.
+        Additional `~imas.db_entry.DBEntry` options, such as ``dd_version`` or ``xml_path``.
 
     Returns
     -------
@@ -116,14 +118,14 @@ def load_bolometers(
 
     If you have a local IMAS database and store the "bolometer.h5" file there:
 
-    >>> bolometers = load_bolometers("imas:hdf5?path=path/to/db/", "r", parent=world)
+    >>> bolometers = load_bolometers("imas:hdf5?path=path/to/db/", parent=world)
 
     If you want to load netCDF files directly:
 
-    >>> bolometers = load_bolometers("path/to/bolometer_file.nc", "r", parent=world)
+    >>> bolometers = load_bolometers("path/to/bolometer_file.nc", parent=world)
     """
     # Load bolometer IDS
-    with DBEntry(*args, **kwargs) as entry:
+    with _open_dbentry_for_reading(*args, **kwargs) as entry:
         # Get available time slices
         ids = get_ids_time_slice(entry, "bolometer")
 
@@ -502,7 +504,7 @@ def visualize(
 
     Examples
     --------
-    >>> bolometers = load_bolometers("imas:hdf5?path=path/to/db/", "r")
+    >>> bolometers = load_bolometers("imas:hdf5?path=path/to/db/")
     >>> fig = visualize(bolometers[0], num_rays=100, ray_from_channel=[0, 3])
     """
     try:
diff --git a/src/cherab/imas/plasma/blend.py b/src/cherab/imas/plasma/blend.py
index 7a3abb0..fc2b781 100644
--- a/src/cherab/imas/plasma/blend.py
+++ b/src/cherab/imas/plasma/blend.py
@@ -31,9 +31,9 @@
 from cherab.core import AtomicData, Maxwellian, Plasma, Species
 from cherab.core.math import VectorAxisymmetricMapper
 from cherab.tools.equilibrium import EFITEquilibrium
-from imas import DBEntry
 from imas.ids_structure import IDSStructure
 
+from .._dbentry import _open_dbentry_for_reading
 from ..ids.common import get_ids_time_slice
 from ..ids.common.ggd import load_grid
 from ..ids.common.grid_radial import get_psi_norm, load_core_grid
@@ -93,17 +93,20 @@ def load_plasma(
     Parameters
     ----------
     *args
-        Arguments passed to the `~imas.db_entry.DBEntry` constructor.
+        IMAS URI, netCDF path, or legacy positional arguments for the core plasma's
+        `~imas.db_entry.DBEntry`. For a URI or path, read mode is selected automatically;
+        do not pass ``"r"``.
     time
         Time for the core plasma, by default 0.
     occurrence_core
         Occurrence index of the ``core_profiles`` IDS, by default 0.
     edge_args
-        Arguments passed to the `~imas.db_entry.DBEntry` constructor for the edge plasma
-        if different from the core plasma. By default None: uses the same as `*args`.
+        URI, netCDF path, or legacy positional DBEntry arguments for the edge plasma if
+        different from the core plasma. Read mode is selected automatically. By default
+        None: uses the same source as `*args`.
     edge_kwargs
-        Keyword arguments passed to the `~imas.db_entry.DBEntry` constructor for the edge plasma
-        if different from the core plasma. By default None: uses the same as `**kwargs`.
+        Additional DBEntry options for the edge plasma if different from the core plasma.
+        By default None: uses the same options as `**kwargs`.
     time_edge
         Time for the edge plasma. If None, uses `~cherab.imas.plasma.load_plasma.time`.
         By default None.
@@ -142,7 +145,8 @@ def load_plasma(
         Parent node in the Raysect scene graph, by default None.
         Typically a `~raysect.optical.scenegraph.world.World` instance.
     **kwargs
-        Keyword arguments passed to the `~imas.db_entry.DBEntry` constructor.
+        Additional `~imas.db_entry.DBEntry` options for the core plasma, such as
+        ``dd_version`` or ``xml_path``.
 
     Returns
     -------
@@ -173,7 +177,7 @@ def load_plasma(
 
     # === Core profiles IDS ===
     try:
-        with DBEntry(*args, **kwargs) as entry_core:
+        with _open_dbentry_for_reading(*args, **kwargs) as entry_core:
             core_profiles_ids = get_ids_time_slice(
                 entry_core,
                 "core_profiles",
@@ -199,7 +203,7 @@ def load_plasma(
 
     # === Edge profiles IDS ===
     try:
-        with DBEntry(*edge_args, **edge_kwargs) as entry_edge:
+        with _open_dbentry_for_reading(*edge_args, **edge_kwargs) as entry_edge:
             edge_profiles_ids = get_ids_time_slice(
                 entry_edge,
                 "edge_profiles",
@@ -277,7 +281,7 @@ def load_plasma(
         not len(grid_ggd.space) and len(grid_ggd.path) and "#" in str(grid_ggd.path)
     )
     if needs_external_grid_reference:
-        with DBEntry(*edge_args, **edge_kwargs) as entry:
+        with _open_dbentry_for_reading(*edge_args, **edge_kwargs) as entry:
             grid, subsets, subset_id = load_grid(grid_ggd, with_subsets=True, entry=entry)
     else:
         grid, subsets, subset_id = load_grid(grid_ggd, with_subsets=True)
diff --git a/src/cherab/imas/plasma/core.py b/src/cherab/imas/plasma/core.py
index c3c7db2..fbca934 100644
--- a/src/cherab/imas/plasma/core.py
+++ b/src/cherab/imas/plasma/core.py
@@ -33,8 +33,8 @@
 from cherab.core.math import VectorAxisymmetricMapper
 from cherab.imas.ids.core_profiles.load_profiles import ProfileData
 from cherab.tools.equilibrium import EFITEquilibrium
-from imas import DBEntry
 
+from .._dbentry import _open_dbentry_for_reading
 from ..ids.common import get_ids_time_slice
 from ..ids.common.grid_radial import get_psi_norm, load_core_grid
 from ..ids.core_profiles import load_core_species
@@ -79,7 +79,9 @@ def load_core_plasma(
     Parameters
     ----------
     *args
-        Arguments passed to the `~imas.db_entry.DBEntry` constructor.
+        IMAS URI, netCDF path, or legacy positional arguments for
+        `~imas.db_entry.DBEntry`. For a URI or path, read mode is selected automatically;
+        do not pass ``"r"``.
     time
         Time for the core plasma, by default 0.
     occurrence
@@ -108,7 +110,7 @@ def load_core_plasma(
         Parent node in the Raysect scene graph, by default None.
         Typically a `~raysect.optical.scenegraph.world.World` instance.
     **kwargs
-        Keyword arguments passed to the `~imas.db_entry.DBEntry` constructor.
+        Additional `~imas.db_entry.DBEntry` options, such as ``dd_version`` or ``xml_path``.
 
     Returns
     -------
@@ -128,7 +130,7 @@ def load_core_plasma(
     # ---------------------------------------------------------
     # Load required data from the core_profiles IDS and form the core grid and species composition
     # data structures.
-    with DBEntry(*args, **kwargs) as entry:
+    with _open_dbentry_for_reading(*args, **kwargs) as entry:
         core_profiles_ids = get_ids_time_slice(
             entry, "core_profiles", time=time, occurrence=occurrence, time_threshold=time_threshold
         )
diff --git a/src/cherab/imas/plasma/edge.py b/src/cherab/imas/plasma/edge.py
index 97d7603..bfe399f 100644
--- a/src/cherab/imas/plasma/edge.py
+++ b/src/cherab/imas/plasma/edge.py
@@ -33,9 +33,9 @@
 from cherab.core import AtomicData, Maxwellian, Plasma, Species
 from cherab.core.math import AxisymmetricMapper, VectorAxisymmetricMapper
 from cherab.tools.equilibrium.efit import FluxSurfaceNormal, PoloidalFieldVector
-from imas import DBEntry
 from imas.ids_structure import IDSStructure
 
+from .._dbentry import _open_dbentry_for_reading
 from ..ggd.base_mesh import GGDGrid
 from ..ids.common import get_ids_time_slice
 from ..ids.common.ggd import load_grid
@@ -89,7 +89,9 @@ def load_edge_plasma(
     Parameters
     ----------
     *args
-        Arguments passed to the `~imas.db_entry.DBEntry` constructor.
+        IMAS URI, netCDF path, or legacy positional arguments for
+        `~imas.db_entry.DBEntry`. For a URI or path, read mode is selected automatically;
+        do not pass ``"r"``.
     time
         Time for the edge plasma, by default 0.
     occurrence
@@ -114,7 +116,7 @@ def load_edge_plasma(
         Parent node in the Raysect scene graph, by default None.
         Typically a `~raysect.optical.scenegraph.world.World` instance.
     **kwargs
-        Keyword arguments passed to the `~imas.db_entry.DBEntry` constructor.
+        Additional `~imas.db_entry.DBEntry` options, such as ``dd_version`` or ``xml_path``.
 
     Returns
     -------
@@ -131,7 +133,7 @@ def load_edge_plasma(
     # -----------------------------------------
     # Load required data from the edge_profiles IDS and form the edge grid and species composition
     # data structures.
-    with DBEntry(*args, **kwargs) as entry:
+    with _open_dbentry_for_reading(*args, **kwargs) as entry:
         edge_profiles_ids = get_ids_time_slice(
             entry, "edge_profiles", time=time, occurrence=occurrence, time_threshold=time_threshold
         )
diff --git a/src/cherab/imas/plasma/equilibrium.py b/src/cherab/imas/plasma/equilibrium.py
index 36fe916..0bab1b3 100644
--- a/src/cherab/imas/plasma/equilibrium.py
+++ b/src/cherab/imas/plasma/equilibrium.py
@@ -26,8 +26,8 @@
 
 from cherab.imas.ids.equilibrium.load_equilibrium import Equilibrium2DData
 from cherab.tools.equilibrium import EFITEquilibrium
-from imas import DBEntry
 
+from .._dbentry import _open_dbentry_for_reading
 from ..ids.common import get_ids_time_slice
 from ..ids.equilibrium import load_equilibrium_data, load_magnetic_field_data
 
@@ -69,7 +69,9 @@ def load_equilibrium(
     Parameters
     ----------
     *args
-        Arguments passed to `~imas.db_entry.DBEntry`.
+        IMAS URI, netCDF path, or legacy positional arguments for
+        `~imas.db_entry.DBEntry`. For a URI or path, read mode is selected automatically;
+        do not pass ``"r"``.
     time
         Time for the equilibrium, by default 0.
     occurrence
@@ -81,7 +83,7 @@ def load_equilibrium(
         If True, returns the ``psi_norm(rho_tor_norm)`` interpolator; otherwise, returns only the
         equilibrium object.
     **kwargs
-        Keyword arguments passed to `~imas.db_entry.DBEntry`.
+        Additional `~imas.db_entry.DBEntry` options, such as ``dd_version`` or ``xml_path``.
 
     Returns
     -------
@@ -98,7 +100,7 @@ def load_equilibrium(
         If the equilibrium IDS does not have a time slice or if ``rho_tor_norm`` is not available
         when ``with_psi_interpolator`` is True.
     """
-    with DBEntry(*args, **kwargs) as entry:
+    with _open_dbentry_for_reading(*args, **kwargs) as entry:
         equilibrium_ids = get_ids_time_slice(
             entry, "equilibrium", time=time, occurrence=occurrence, time_threshold=time_threshold
         )
@@ -159,7 +161,9 @@ def load_magnetic_field(
     Parameters
     ----------
     *args
-        Arguments passed to `~imas.db_entry.DBEntry`.
+        IMAS URI, netCDF path, or legacy positional arguments for
+        `~imas.db_entry.DBEntry`. For a URI or path, read mode is selected automatically;
+        do not pass ``"r"``.
     time
         Time for the equilibrium, by default 0.
     occurrence
@@ -168,7 +172,7 @@ def load_magnetic_field(
         Maximum allowed difference between the requested time and the nearest available
         time, by default `numpy.inf`.
     **kwargs
-        Keyword arguments passed to `~imas.db_entry.DBEntry`.
+        Additional `~imas.db_entry.DBEntry` options, such as ``dd_version`` or ``xml_path``.
 
     Returns
     -------
@@ -180,7 +184,7 @@ def load_magnetic_field(
     RuntimeError
         If the equilibrium IDS does not have a time slice.
     """
-    with DBEntry(*args, **kwargs) as entry:
+    with _open_dbentry_for_reading(*args, **kwargs) as entry:
         equilibrium_ids = get_ids_time_slice(
             entry, "equilibrium", time=time, occurrence=occurrence, time_threshold=time_threshold
         )
diff --git a/src/cherab/imas/wall/wall.py b/src/cherab/imas/wall/wall.py
index 4ea2956..06ec4de 100644
--- a/src/cherab/imas/wall/wall.py
+++ b/src/cherab/imas/wall/wall.py
@@ -23,8 +23,7 @@
 from raysect.optical.material.material import Material
 from raysect.primitive import Mesh
 
-from imas import DBEntry
-
+from .._dbentry import _open_dbentry_for_reading
 from ..ids.common import get_ids_time_slice
 from ..ids.wall import load_wall_2d, load_wall_3d
 
@@ -47,7 +46,9 @@ def load_wall_mesh(
     Parameters
     ----------
     *args
-        Arguments passed to the `~imas.db_entry.DBEntry` constructor.
+        IMAS URI, netCDF path, or legacy positional arguments for
+        `~imas.db_entry.DBEntry`. For a URI or path, read mode is selected automatically;
+        do not pass ``"r"``.
     time
         Time for the wall, by default 0.
     occurrence
@@ -70,7 +71,7 @@ def load_wall_mesh(
         Parent node in the Raysect scene graph, by default None.
         Normally, `~raysect.optical.scenegraph.world.World` instance.
     **kwargs
-        Keyword arguments passed to the `~imas.db_entry.DBEntry` constructor.
+        Additional `~imas.db_entry.DBEntry` options, such as ``dd_version`` or ``xml_path``.
 
     Returns
     -------
@@ -82,12 +83,12 @@ def load_wall_mesh(
     >>> from raysect.optical import World
     >>> world = World()
     >>> meshes = load_wall_mesh(
-    ...     "imas:hdf5?path=/work/imas/shared/imasdb/ITER_MD/3/116100/1001/", "r", parent=world
+    ...     "imas:hdf5?path=/work/imas/shared/imasdb/ITER_MD/3/116100/1001/", parent=world
     ... )
     >>> meshes
     {'FullTokamak.none.none': <raysect.primitive.mesh.mesh.Mesh at 0x1766322a0>}
     """
-    with DBEntry(*args, **kwargs) as entry:
+    with _open_dbentry_for_reading(*args, **kwargs) as entry:
         wall_ids = get_ids_time_slice(
             entry, "wall", time=time, occurrence=occurrence, time_threshold=time_threshold
         )
@@ -115,13 +116,15 @@ def load_wall_outline(
     Parameters
     ----------
     *args
-        Arguments passed to the `~imas.db_entry.DBEntry` constructor.
+        IMAS URI, netCDF path, or legacy positional arguments for
+        `~imas.db_entry.DBEntry`. For a URI or path, read mode is selected automatically;
+        do not pass ``"r"``.
     occurrence
         Occurrence index of the ``wall`` IDS, by default 0.
     desc_index
         Index of ``description_2d``, by default 0.
     **kwargs
-        Keyword arguments passed to the `~imas.db_entry.DBEntry` constructor.
+        Additional `~imas.db_entry.DBEntry` options, such as ``dd_version`` or ``xml_path``.
 
     Returns
     -------
@@ -130,7 +133,7 @@ def load_wall_outline(
 
     Examples
     --------
-    >>> load_wall_outline("imas:hdf5?path=/work/imas/shared/imasdb/ITER_MD/3/116000/5/", "r")
+    >>> load_wall_outline("imas:hdf5?path=/work/imas/shared/imasdb/ITER_MD/3/116000/5/")
     {'First Wall': array([[ 4.11129713, -2.49559808],
                           [ 4.11129713, -1.48329401],
                           ...
@@ -139,7 +142,7 @@ def load_wall_outline(
                         ...
                         [ 6.36320019, -3.24460006]])}
     """
-    with DBEntry(*args, **kwargs) as entry:
+    with _open_dbentry_for_reading(*args, **kwargs) as entry:
         description2d = entry.get("wall", occurrence=occurrence, autoconvert=False).description_2d[
             desc_index
         ]

From 25f5ab69ed02a47934922bdc666ea71b2676ce6e Mon Sep 17 00:00:00 2001
From: munechika-koyo <munechika.koyo@gmail.com>
Date: Tue, 11 Aug 2026 11:50:53 +0200
Subject: [PATCH 24/36] =?UTF-8?q?=F0=9F=94=A5=20Remove=20"r"=20mode=20para?=
 =?UTF-8?q?meter?=
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit

---
 demos/ITER/core_plasma_plot_profiles.py     |  4 +-
 demos/ITER/edge_plasma_plot_profiles.py     | 20 ++++--
 demos/ITER/full_plasma_plot_profiles.py     | 18 ++++--
 docs/notebooks/observer/bolometer.ipynb     |  2 +-
 docs/notebooks/plasma/1_equilibrium.ipynb   |  2 +-
 docs/notebooks/plasma/3_full_plasma.ipynb   |  2 +-
 docs/notebooks/plasma/4_emission.ipynb      |  2 +-
 docs/notebooks/radiation/radiation_3d.ipynb |  1 -
 docs/source/intro.md                        | 27 +++++++-
 tests/emitter/test_radiation.py             |  6 +-
 tests/observer/test_bolometer.py            |  2 +-
 tests/plasma/test_blend.py                  |  2 +-
 tests/plasma/test_core.py                   |  2 +-
 tests/plasma/test_edge.py                   |  2 +-
 tests/plasma/test_equilibrium.py            |  4 +-
 tests/test_dbentry.py                       | 71 +++++++++++++++++++++
 16 files changed, 139 insertions(+), 28 deletions(-)
 create mode 100644 tests/test_dbentry.py

diff --git a/demos/ITER/core_plasma_plot_profiles.py b/demos/ITER/core_plasma_plot_profiles.py
index e4a05ea..105e5cd 100644
--- a/demos/ITER/core_plasma_plot_profiles.py
+++ b/demos/ITER/core_plasma_plot_profiles.py
@@ -73,7 +73,7 @@ def plot_quantity(quantity, extent, title="", logscale=False, symmetric=False) -
 
 # Load and plot equilibrium
 path = iter_jintrac()
-equilibrium, psi_interpolator = load_equilibrium(path, "r", with_psi_interpolator=True)
+equilibrium, psi_interpolator = load_equilibrium(path, with_psi_interpolator=True)
 plot_equilibrium(equilibrium)
 plt.gcf().savefig(plots_path / "equilibrium.png", dpi=200)
 
@@ -95,7 +95,7 @@ def plot_quantity(quantity, extent, title="", logscale=False, symmetric=False) -
 fig.savefig(plots_path / "bz.png", dpi=200)
 
 # Load core plasma
-plasma = load_core_plasma(path, "r", equilibrium=equilibrium, psi_interpolator=psi_interpolator)
+plasma = load_core_plasma(path, equilibrium=equilibrium, psi_interpolator=psi_interpolator)
 
 # Sample and plot electron profiles
 ne_plasma = sample3d_grid(plasma.electron_distribution.density, xsamp, [0], zsamp)
diff --git a/demos/ITER/edge_plasma_plot_profiles.py b/demos/ITER/edge_plasma_plot_profiles.py
index 1dda2ab..2eeba93 100644
--- a/demos/ITER/edge_plasma_plot_profiles.py
+++ b/demos/ITER/edge_plasma_plot_profiles.py
@@ -79,18 +79,22 @@ def plot_quantity(quantity, extent, title="", logscale=False, symmetric=False) -
 
 # Load and plot equilibrium
 try:
-    equilibrium = load_equilibrium(iter_jintrac(), "r")
+    equilibrium = load_equilibrium(iter_jintrac())
     plot_equilibrium(equilibrium)
     plt.gcf().savefig(plots_path / "equilibrium.png", dpi=200)
     b_field = equilibrium.b_field
 except RuntimeError:
     try:
-        b_field = load_magnetic_field(iter_jintrac(), "r")
+        b_field = load_magnetic_field(iter_jintrac())
     except RuntimeError:
         b_field = None
 
 # Sample and plot magnetic field
 plot_velocity = False
+b = None
+b_length = None
+radial_vector = None
+poloidal_vector = None
 if b_field is not None:
     try:
         xsamp, zsamp, b = samplevector2d(b_field, (xl, xu, nx), (zl, zu, nz))
@@ -116,7 +120,7 @@ def plot_quantity(quantity, extent, title="", logscale=False, symmetric=False) -
         )
 
 # Load edge plasma
-plasma = load_edge_plasma(iter_jintrac(), "r", b_field=b_field)
+plasma = load_edge_plasma(iter_jintrac(), b_field=b_field)
 
 # Sample and plot electron profiles
 xsamp, _, zsamp, ne_plasma = sample3d(
@@ -132,7 +136,7 @@ def plot_quantity(quantity, extent, title="", logscale=False, symmetric=False) -
 fig = plot_quantity(te_plasma, extent, title="Te [eV]", logscale=True)
 fig.savefig(plots_path / "edge_te.png", dpi=200)
 
-if plot_velocity:
+if plot_velocity and b is not None and b_length is not None:
     electron_velocity = samplevector3d_grid(
         plasma.electron_distribution.bulk_velocity, xsamp, [0], zsamp
     ).squeeze()
@@ -174,7 +178,13 @@ def plot_quantity(quantity, extent, title="", logscale=False, symmetric=False) -
             dpi=200,
         )
 
-        if plot_velocity:
+        if (
+            plot_velocity
+            and b is not None
+            and b_length is not None
+            and radial_vector is not None
+            and poloidal_vector is not None
+        ):
             velocity = samplevector3d_grid(
                 species.distribution.bulk_velocity, xsamp, [0], zsamp
             ).squeeze()
diff --git a/demos/ITER/full_plasma_plot_profiles.py b/demos/ITER/full_plasma_plot_profiles.py
index 6e802ce..6ae1c11 100644
--- a/demos/ITER/full_plasma_plot_profiles.py
+++ b/demos/ITER/full_plasma_plot_profiles.py
@@ -79,13 +79,17 @@ def plot_quantity(quantity, extent, title="", logscale=False, symmetric=False) -
 extent = [xl, xu, zl, zu]
 
 # Load and plot equilibrium
-equilibrium, psi_interpolator = load_equilibrium(iter_jintrac(), "r", with_psi_interpolator=True)
+equilibrium, psi_interpolator = load_equilibrium(iter_jintrac(), with_psi_interpolator=True)
 plot_equilibrium(equilibrium)
 plt.gcf().savefig(plots_path / "equilibrium.png", dpi=200)
 
 # Sample and plot magnetic field
 plot_velocity = False
 b_field = equilibrium.b_field
+b = None
+b_length = None
+radial_vector = None
+poloidal_vector = None
 try:
     xsamp, zsamp, b = samplevector2d(b_field, (xl, xu, nx), (zl, zu, nz))
     b_length = np.sqrt((b * b).sum(2))
@@ -110,7 +114,7 @@ def plot_quantity(quantity, extent, title="", logscale=False, symmetric=False) -
     )
 
 # Load edge plasma
-plasma = load_plasma(iter_jintrac(), "r", equilibrium=equilibrium)
+plasma = load_plasma(iter_jintrac(), equilibrium=equilibrium)
 
 # Sample and plot electron profiles
 xsamp, _, zsamp, ne_plasma = sample3d(
@@ -126,7 +130,7 @@ def plot_quantity(quantity, extent, title="", logscale=False, symmetric=False) -
 fig = plot_quantity(te_plasma, extent, title="Te [eV]", logscale=True)
 fig.savefig(plots_path / "te.png", dpi=200)
 
-if plot_velocity:
+if plot_velocity and b is not None and b_length is not None:
     electron_velocity = samplevector3d_grid(
         plasma.electron_distribution.bulk_velocity, xsamp, [0], zsamp
     ).squeeze()
@@ -168,7 +172,13 @@ def plot_quantity(quantity, extent, title="", logscale=False, symmetric=False) -
             dpi=200,
         )
 
-        if plot_velocity:
+        if (
+            plot_velocity
+            and b is not None
+            and b_length is not None
+            and radial_vector is not None
+            and poloidal_vector is not None
+        ):
             velocity = samplevector3d_grid(
                 species.distribution.bulk_velocity, xsamp, [0], zsamp
             ).squeeze()
diff --git a/docs/notebooks/observer/bolometer.ipynb b/docs/notebooks/observer/bolometer.ipynb
index dd73352..e1385c4 100644
--- a/docs/notebooks/observer/bolometer.ipynb
+++ b/docs/notebooks/observer/bolometer.ipynb
@@ -52,7 +52,7 @@
    "outputs": [],
    "source": [
     "path = bolometer_moc()\n",
-    "bolometers = load_bolometers(path, \"r\")"
+    "bolometers = load_bolometers(path)"
    ]
   },
   {
diff --git a/docs/notebooks/plasma/1_equilibrium.ipynb b/docs/notebooks/plasma/1_equilibrium.ipynb
index 75ad162..63148eb 100644
--- a/docs/notebooks/plasma/1_equilibrium.ipynb
+++ b/docs/notebooks/plasma/1_equilibrium.ipynb
@@ -59,7 +59,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "equilibrium = load_equilibrium(path, \"r\")"
+    "equilibrium = load_equilibrium(path)"
    ]
   },
   {
diff --git a/docs/notebooks/plasma/3_full_plasma.ipynb b/docs/notebooks/plasma/3_full_plasma.ipynb
index 2fa643b..68e56c5 100644
--- a/docs/notebooks/plasma/3_full_plasma.ipynb
+++ b/docs/notebooks/plasma/3_full_plasma.ipynb
@@ -203,7 +203,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "plasma = load_plasma(path, \"r\")"
+    "plasma = load_plasma(path)"
    ]
   },
   {
diff --git a/docs/notebooks/plasma/4_emission.ipynb b/docs/notebooks/plasma/4_emission.ipynb
index 4c89862..3b01fe7 100644
--- a/docs/notebooks/plasma/4_emission.ipynb
+++ b/docs/notebooks/plasma/4_emission.ipynb
@@ -127,7 +127,7 @@
    "source": [
     "world = World()\n",
     "path = iter_jintrac()\n",
-    "plasma = load_plasma(path, \"r\", parent=world, atomic_data=atomic_data)"
+    "plasma = load_plasma(path, parent=world, atomic_data=atomic_data)"
    ]
   },
   {
diff --git a/docs/notebooks/radiation/radiation_3d.ipynb b/docs/notebooks/radiation/radiation_3d.ipynb
index 9b66f18..e79e36e 100644
--- a/docs/notebooks/radiation/radiation_3d.ipynb
+++ b/docs/notebooks/radiation/radiation_3d.ipynb
@@ -68,7 +68,6 @@
     "world = World()\n",
     "emitter = load_radiation_emitter(\n",
     "    path,\n",
-    "    \"r\",\n",
     "    time=0.0042853,\n",
     "    parent=world,\n",
     "    interpolator_cache=\"disk\",\n",
diff --git a/docs/source/intro.md b/docs/source/intro.md
index cc7377d..6420595 100644
--- a/docs/source/intro.md
+++ b/docs/source/intro.md
@@ -60,16 +60,37 @@ from cherab.imas.plasma import load_plasma
 world = World()
 
 # Load plasma from IMAS database
-# You can put same parameters defined in `imas.DBEntry`
 plasma = load_plasma(
     "imas:hdf5?path=testdb",  # IMAS URI
-    "r",  # read mode
     time=0.0,
     parent=world,
 )
 ```
 
-This script creates a `World` object and loads a plasma from an IMAS database located at `testdb`. The plasma is created close to time `0.0` and is added to the world as its parent.
+This script creates a `World` object and loads a plasma from an IMAS database located at `testdb`, reading the `core_profiles/equilibrium/edge_profiles` IDSs.
+The plasma is created close to time `0.0` and is added to the world as its parent.
+
+All `cherab-imas` loader APIs open URI and netCDF sources in IMAS read mode automatically. Pass
+the source directly, without the `"r"` required by `imas.DBEntry` itself:
+
+```python
+# IMAS URI
+plasma = load_plasma("imas:hdf5?path=testdb")
+
+# A netCDF data entry works the same way
+plasma = load_plasma("/path/to/data.nc")
+
+# Other DBEntry options remain available as keyword arguments
+plasma = load_plasma("/path/to/data.nc", dd_version="4.1.0")
+
+# Legacy form is also supported
+from imas.ids_defs import HDF5_BACKEND
+
+plasma = load_plasma(HDF5_BACKEND, "testdb", 12345, 0)
+```
+
+Existing calls that explicitly pass `"r"` continue to work for compatibility, but are deprecated.
+Modes that create or replace data, such as `"w"`, `"a"`, and `"x"`, are rejected by loader APIs.
 
 You can find more examples and detailed documentation in the [Examples](examples) and [API Reference](api) sections.
 
diff --git a/tests/emitter/test_radiation.py b/tests/emitter/test_radiation.py
index 29da939..c01da03 100644
--- a/tests/emitter/test_radiation.py
+++ b/tests/emitter/test_radiation.py
@@ -10,6 +10,7 @@
 from raysect.primitive import Cylinder, Subtract
 
 import cherab.imas.emitter.radiation as radiation_module
+from cherab.imas import _dbentry as dbentry_module
 from cherab.imas.emitter import load_radiation_emitter
 from cherab.imas.plasma.equilibrium import load_equilibrium
 
@@ -63,14 +64,14 @@ def _patch_dbentry_for_open_entries(
     monkeypatch: pytest.MonkeyPatch,
     entries: dict[tuple, DBEntry],
 ) -> None:
-    original_dbentry = radiation_module.DBEntry
+    original_dbentry = dbentry_module.DBEntry
 
     def _dbentry_router(*args, **kwargs):
         if not kwargs and args in entries:
             return _OpenEntryContext(entries[args])
         return original_dbentry(*args, **kwargs)
 
-    monkeypatch.setattr(radiation_module, "DBEntry", _dbentry_router)
+    monkeypatch.setattr(dbentry_module, "DBEntry", _dbentry_router)
 
 
 def _write_split_radiation_to_memory(
@@ -146,7 +147,6 @@ def average_gaussian_faces_per_toroidal(self, phis):
 
     primitive = load_radiation_emitter(
         path_iter_jorek,
-        "r",
         source="coefficients",
         **_cache_kwargs(radiation_interpolator_cache),
     )
diff --git a/tests/observer/test_bolometer.py b/tests/observer/test_bolometer.py
index 244e2aa..6d704a6 100644
--- a/tests/observer/test_bolometer.py
+++ b/tests/observer/test_bolometer.py
@@ -10,7 +10,7 @@
 
 def test_load_bolometers(path_bolometer_moc: str) -> None:
     """Test loading bolometer data from an IDS dataset."""
-    bolometers = load_bolometers(path_bolometer_moc, "r")
+    bolometers = load_bolometers(path_bolometer_moc)
 
     # Check that the bolometer cameras are loaded correctly
     assert len(bolometers) == 3
diff --git a/tests/plasma/test_blend.py b/tests/plasma/test_blend.py
index 774b629..44a691a 100644
--- a/tests/plasma/test_blend.py
+++ b/tests/plasma/test_blend.py
@@ -10,7 +10,7 @@
 
 def test_load_plasma(path_iter_jintrac: str):
     """Test basic loading of plasma data from an IMAS file."""
-    plasma = load_plasma(path_iter_jintrac, "r")
+    plasma = load_plasma(path_iter_jintrac)
 
     # Test that a Plasma object is returned
     assert isinstance(plasma, Plasma)
diff --git a/tests/plasma/test_core.py b/tests/plasma/test_core.py
index ff23edb..4bbbc4d 100644
--- a/tests/plasma/test_core.py
+++ b/tests/plasma/test_core.py
@@ -10,7 +10,7 @@
 
 def test_load_core_plasma(path_iter_jintrac: str):
     """Test basic loading of core plasma data from an IMAS file."""
-    plasma = load_core_plasma(path_iter_jintrac, "r")
+    plasma = load_core_plasma(path_iter_jintrac)
 
     # Test that a Plasma object is returned
     assert isinstance(plasma, Plasma)
diff --git a/tests/plasma/test_edge.py b/tests/plasma/test_edge.py
index 51998be..8009cd7 100644
--- a/tests/plasma/test_edge.py
+++ b/tests/plasma/test_edge.py
@@ -19,7 +19,7 @@
 
 def test_load_edge_plasma(path_iter_jintrac: str):
     """Test loading of edge plasma data from an IMAS file."""
-    plasma = load_edge_plasma(path_iter_jintrac, "r")
+    plasma = load_edge_plasma(path_iter_jintrac)
 
     # Test that a Plasma object is returned
     assert isinstance(plasma, Plasma)
diff --git a/tests/plasma/test_equilibrium.py b/tests/plasma/test_equilibrium.py
index ccb1f65..3bf538a 100644
--- a/tests/plasma/test_equilibrium.py
+++ b/tests/plasma/test_equilibrium.py
@@ -9,7 +9,7 @@
 
 def test_load_equilibrium(path_iter_jintrac: str):
     """Test loading of equilibrium data from an IMAS file."""
-    equilibrium = load_equilibrium(path_iter_jintrac, "r")
+    equilibrium = load_equilibrium(path_iter_jintrac)
 
     # Test that equilibrium object is returned
     assert isinstance(equilibrium, EFITEquilibrium)
@@ -20,7 +20,7 @@ def test_load_equilibrium(path_iter_jintrac: str):
 
 def test_load_magnetic_field(path_iter_jintrac: str):
     """Test loading of magnetic field data from an IMAS file."""
-    magnetic_field = load_magnetic_field(path_iter_jintrac, "r")
+    magnetic_field = load_magnetic_field(path_iter_jintrac)
 
     # Test that Vector3DFunction2D object is returned
     assert isinstance(magnetic_field, Vector3DFunction2D)
diff --git a/tests/test_dbentry.py b/tests/test_dbentry.py
new file mode 100644
index 0000000..b6e3fc5
--- /dev/null
+++ b/tests/test_dbentry.py
@@ -0,0 +1,71 @@
+from pathlib import Path
+
+import pytest
+
+from cherab.imas import _dbentry
+
+
+class _DummyDBEntry:
+    pass
+
+
+@pytest.fixture
+def dbentry_spy(monkeypatch: pytest.MonkeyPatch) -> list[tuple[tuple, dict]]:
+    calls = []
+
+    def create_dbentry(*args, **kwargs):
+        calls.append((args, kwargs))
+        return _DummyDBEntry()
+
+    monkeypatch.setattr(_dbentry, "DBEntry", create_dbentry)
+    return calls
+
+
+@pytest.mark.parametrize("uri", ["imas:hdf5?path=testdb", Path("data.nc")])
+def test_uri_mode_is_added_automatically(uri, dbentry_spy):
+    entry = _dbentry._open_dbentry_for_reading(uri, dd_version="4.1.0")
+
+    assert isinstance(entry, _DummyDBEntry)
+    assert dbentry_spy == [((str(uri), "r"), {"dd_version": "4.1.0"})]
+
+
+def test_keyword_uri_mode_is_added_automatically(dbentry_spy):
+    _dbentry._open_dbentry_for_reading(uri="data.nc", xml_path="IDSDef.xml")
+
+    assert dbentry_spy == [(("data.nc", "r"), {"xml_path": "IDSDef.xml"})]
+
+
+@pytest.mark.parametrize(
+    ("args", "kwargs"),
+    [
+        (("data.nc", "r"), {}),
+        (("data.nc",), {"mode": "r"}),
+        ((), {"uri": "data.nc", "mode": "r"}),
+    ],
+)
+def test_explicit_read_mode_is_deprecated(args, kwargs, dbentry_spy):
+    with pytest.deprecated_call(match="selected automatically"):
+        _dbentry._open_dbentry_for_reading(*args, **kwargs)
+
+    assert dbentry_spy == [(("data.nc", "r"), {})]
+
+
+@pytest.mark.parametrize("mode", ["w", "a", "x", "r+"])
+def test_non_read_mode_is_rejected(mode, dbentry_spy):
+    with pytest.raises(ValueError, match="only support mode 'r'"):
+        _dbentry._open_dbentry_for_reading("data.nc", mode)
+
+    assert not dbentry_spy
+
+
+def test_duplicate_mode_is_rejected(dbentry_spy):
+    with pytest.raises(TypeError, match="both positionally and by keyword"):
+        _dbentry._open_dbentry_for_reading("data.nc", "r", mode="r")
+
+    assert not dbentry_spy
+
+
+def test_legacy_constructor_is_passed_through(dbentry_spy):
+    _dbentry._open_dbentry_for_reading(13, "ITER", 12345, 1, "user", dd_version="3.42.0")
+
+    assert dbentry_spy == [((13, "ITER", 12345, 1, "user"), {"dd_version": "3.42.0"})]

From bab3a7f3486f0488983326bdcdad2770b9b14dc9 Mon Sep 17 00:00:00 2001
From: munechika-koyo <munechika.koyo@gmail.com>
Date: Tue, 11 Aug 2026 11:51:40 +0200
Subject: [PATCH 25/36] =?UTF-8?q?=F0=9F=94=A5=20Remove=20unused=20phis=20p?=
 =?UTF-8?q?arameter?=
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit

---
 src/cherab/imas/emitter/radiation.py | 16 ++++---------
 tests/emitter/test_radiation.py      | 36 +++-------------------------
 2 files changed, 7 insertions(+), 45 deletions(-)

diff --git a/src/cherab/imas/emitter/radiation.py b/src/cherab/imas/emitter/radiation.py
index 49208b4..bd59f73 100644
--- a/src/cherab/imas/emitter/radiation.py
+++ b/src/cherab/imas/emitter/radiation.py
@@ -211,7 +211,6 @@ def load_radiation_emitter(
     psi_interpolator: Callable[[float], float] | None = None,
     mask: Function2D | Function3D | None = None,
     num_toroidal: int = 64,
-    phis: NDArray[np.float64] | None = None,
     source: Literal["auto", "values", "coefficients"] = "auto",
     time_threshold: float = np.inf,
     step: float = 0.01,
@@ -266,8 +265,7 @@ def load_radiation_emitter(
         Radiation process identifier index (or indices) to load.
         By default, all available processes are summed together.
         Reference: https://imas-data-dictionary.readthedocs.io/en/latest/generated/identifier/radiation_identifier.html
-        .. note::
-            The emissivity value array is assumed to follow the same x-axis as the grid subset.
+        The emissivity value array is assumed to follow the same x-axis as the grid subset.
     grid_ggd
         Specific grid GGD structure alternative to the one in the IDS.
     grid_subset_id
@@ -287,9 +285,6 @@ def load_radiation_emitter(
     num_toroidal
         Number of toroidal subdivisions for 3D grid extension, by default 64.
         This is used only when the grid is loaded by `.load_unstruct_grid_2d_extended`.
-    phis
-        Array of toroidal angles in degrees for emissivity reconstruction, by default None.
-        This is used only when the grid is loaded by `.load_unstruct_grid_2d_extended`.
     source
         Source for emissivity data: ``"auto"`` (tries values then coefficients), ``"values"``
         (blended emissivity from core profiles + (edge) GGD values), or ``"coefficients"``
@@ -553,13 +548,10 @@ def load_radiation_emitter(
 
         constructor = FourierBezierConstructor(grid_ggd, coefficients=coeff)
 
-        if phis is None:
-            d_phi = 360.0 / grid.num_toroidal
-            phis_array = np.arange(d_phi * 0.5, 360.0, d_phi, dtype=np.float64)
-        else:
-            phis_array = np.asarray(phis, dtype=np.float64)
+        d_phi = 360.0 / grid.num_toroidal
+        phis = np.arange(d_phi * 0.5, 360.0, d_phi, dtype=np.float64)
 
-        emissivity = constructor.average_gaussian_faces_per_toroidal(phis_array).ravel()
+        emissivity = constructor.average_gaussian_faces_per_toroidal(phis).ravel()
         primitive_name = f"RadiationEmitter_{ids.time[0]}s, entry {entry_reference}"
 
         rad_func = grid.interpolator(
diff --git a/tests/emitter/test_radiation.py b/tests/emitter/test_radiation.py
index c01da03..259c3b0 100644
--- a/tests/emitter/test_radiation.py
+++ b/tests/emitter/test_radiation.py
@@ -161,6 +161,8 @@ def average_gaussian_faces_per_toroidal(self, phis):
     assert np.all(np.diff(phis_used) > 0)
     assert 0.0 < phis_used[0] < 360.0
     assert 0.0 < phis_used[-1] < 360.0
+    d_phi = 360.0 / phis_used.size
+    np.testing.assert_allclose(phis_used, np.arange(0.5 * d_phi, 360.0, d_phi, dtype=np.float64))
 
 
 def test_load_radiation_emitter_auto_falls_back_to_coefficients(
@@ -178,38 +180,6 @@ def test_load_radiation_emitter_auto_falls_back_to_coefficients(
     assert primitive.material is not None
 
 
-def test_load_radiation_emitter_coefficients_uses_given_phis(
-    path_iter_jorek: str,
-    monkeypatch: pytest.MonkeyPatch,
-    radiation_interpolator_cache: tuple[Literal["memory", "disk"], Path | None],
-):
-    captured: dict[str, np.ndarray] = {}
-    original_constructor = radiation_module.FourierBezierConstructor
-
-    class _SpyFourierBezierConstructor:
-        def __init__(self, *args, **kwargs):
-            self._inner = original_constructor(*args, **kwargs)
-
-        def average_gaussian_faces_per_toroidal(self, phis):
-            captured["phis"] = np.asarray(phis, dtype=np.float64).copy()
-            return self._inner.average_gaussian_faces_per_toroidal(phis)
-
-    monkeypatch.setattr(radiation_module, "FourierBezierConstructor", _SpyFourierBezierConstructor)
-
-    expected_phis = np.array([15.0, 105.0, 195.0, 285.0], dtype=np.float64)
-    primitive = load_radiation_emitter(
-        path_iter_jorek,
-        "r",
-        source="coefficients",
-        phis=expected_phis,
-        **_cache_kwargs(radiation_interpolator_cache),
-    )
-
-    assert isinstance(primitive, (Subtract, Cylinder))
-    assert primitive.material is not None
-    np.testing.assert_allclose(captured["phis"], expected_phis)
-
-
 def test_load_radiation_emitter_values_raises_for_jorek(
     path_iter_jorek: str,
     radiation_interpolator_cache: tuple[Literal["memory", "disk"], Path | None],
@@ -430,7 +400,7 @@ class _UnexpectedDBEntry:
         def __init__(self, *args, **kwargs):
             raise AssertionError("DBEntry must not be instantiated for invalid source.")
 
-    monkeypatch.setattr(radiation_module, "DBEntry", _UnexpectedDBEntry)
+    monkeypatch.setattr(dbentry_module, "DBEntry", _UnexpectedDBEntry)
 
     with pytest.raises(
         ValueError,

From 6ea5cfd835dc6ec5a2dd9955521eb2875931a1d3 Mon Sep 17 00:00:00 2001
From: munechika-koyo <munechika.koyo@gmail.com>
Date: Tue, 11 Aug 2026 11:52:30 +0200
Subject: [PATCH 26/36] =?UTF-8?q?=F0=9F=93=9D=20Update=20changelog=20inclu?=
 =?UTF-8?q?ding=20"r"=20rejection?=
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit

---
 CHANGELOG.md | 4 ++++
 1 file changed, 4 insertions(+)

diff --git a/CHANGELOG.md b/CHANGELOG.md
index 51a012f..2e79ac1 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -26,6 +26,10 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 
 ### Changed
 
+- Select read mode automatically when CHERAB object loaders receive an IMAS URI or netCDF path;
+  explicit `"r"` arguments remain temporarily supported but are deprecated
+- Reject data-creating DBEntry modes in CHERAB object loader APIs
+- Update API docstrings, examples, demos, and notebooks to use mode-free loader calls
 - Improve radiation emitter loading checks for duplicate emissivity values and core-profile grid data
 - Refactor `load_grid` to support explicit `entry` selection and improve referenced-grid source resolution
 - Enhance plasma and emitter loading paths to use entry-reference aware grid resolution

From d5547e29aff31fab467dd557657a7efd33073dab Mon Sep 17 00:00:00 2001
From: munechika-koyo <munechika.koyo@gmail.com>
Date: Tue, 11 Aug 2026 11:56:51 +0200
Subject: [PATCH 27/36] =?UTF-8?q?=F0=9F=93=9D=20Add=20radiation=202d=20not?=
 =?UTF-8?q?ebook?=
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit

---
 docs/notebooks/radiation/radiation_2d.ipynb | 192 ++++++++++++++++++++
 1 file changed, 192 insertions(+)
 create mode 100644 docs/notebooks/radiation/radiation_2d.ipynb

diff --git a/docs/notebooks/radiation/radiation_2d.ipynb b/docs/notebooks/radiation/radiation_2d.ipynb
new file mode 100644
index 0000000..5d7600b
--- /dev/null
+++ b/docs/notebooks/radiation/radiation_2d.ipynb
@@ -0,0 +1,192 @@
+{
+ "cells": [
+  {
+   "cell_type": "markdown",
+   "id": "0",
+   "metadata": {},
+   "source": [
+    "# 2D Radiation Emitter\n",
+    "\n",
+    "This notebook demonstrates how to create a radiation emitter from 2D profiles stored in an IMAS [radiation IDS](https://imas-data-dictionary.readthedocs.io/en/latest/generated/ids/radiation.html).\n",
+    "\n",
+    "The example test data was calculated by SOLPS-ITER for an ITER scenario."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "1",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "import numpy as np\n",
+    "import ultraplot as uplt\n",
+    "from raysect.optical import World\n",
+    "\n",
+    "from cherab.imas.datasets import iter_solps\n",
+    "from cherab.imas.emitter import load_radiation_emitter\n",
+    "\n",
+    "# Set dark background for plots\n",
+    "uplt.rc.style = \"dark_background\""
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "2",
+   "metadata": {},
+   "source": [
+    "## Retrieve SOLPS-ITER sample data"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "3",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "path = iter_solps()"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "4",
+   "metadata": {},
+   "source": [
+    "## Create 3D radiation emitter from IMAS IDS"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "5",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "world = World()\n",
+    "emitter = load_radiation_emitter(\n",
+    "    path,\n",
+    "    parent=world,\n",
+    "    interpolator_cache=\"disk\",\n",
+    ")"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "6",
+   "metadata": {},
+   "source": [
+    "## Visualize the emitter in 2D\n",
+    "\n",
+    "Sample 2D visualization of the radiation function."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "7",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "R_MIN, R_MAX = 4.0, 8.5\n",
+    "Z_MIN, Z_MAX = -4.7, 4.8\n",
+    "RES = 0.01  # resolution of grid in [m]\n",
+    "n_r = round((R_MAX - R_MIN) / RES) + 1\n",
+    "n_z = round((Z_MAX - Z_MIN) / RES) + 1\n",
+    "dr, dz = (R_MAX - R_MIN) / (n_r - 1), (Z_MAX - Z_MIN) / (n_z - 1)\n",
+    "\n",
+    "# (r, z) coordinates at phi=0\n",
+    "r_pts = np.linspace(R_MIN, R_MAX, n_r, endpoint=True)\n",
+    "z_pts = np.linspace(Z_MIN, Z_MAX, n_z, endpoint=True)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "8",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "rad = np.zeros((n_r, n_z), dtype=float)\n",
+    "for i, j in np.ndindex(n_r, n_z):\n",
+    "    rad[i, j] = emitter.material.radiation_function(\n",
+    "        r_pts[i],\n",
+    "        0.0,\n",
+    "        z_pts[j],\n",
+    "    )"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "9",
+   "metadata": {},
+   "source": [
+    "### Cross-sectional view of the radiation\n",
+    "\n",
+    "2D poloidal cross-sections of the radiation function."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "10",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "rad[rad <= 0] = np.nan  # set non-positive values to NaN for log scale plotting\n",
+    "\n",
+    "fig, ax = uplt.subplots()\n",
+    "\n",
+    "im = ax.pcolormesh(\n",
+    "    r_pts,\n",
+    "    z_pts,\n",
+    "    rad.T,\n",
+    "    shading=\"auto\",\n",
+    "    cmap=\"inferno\",\n",
+    "    discrete=False,\n",
+    "    norm=\"log\",\n",
+    ")\n",
+    "ax.colorbar(\n",
+    "    im,\n",
+    "    loc=\"r\",\n",
+    "    label=\"[W/m³]\",\n",
+    "    tickminor=True,\n",
+    "    formatter=\"log\",\n",
+    ")\n",
+    "ax.format(\n",
+    "    aspect=\"equal\",\n",
+    "    titleborder=False,\n",
+    "    xlim=(R_MIN, R_MAX),\n",
+    "    ylim=(Z_MIN, Z_MAX),\n",
+    "    xlabel=\"$R$ [m]\",\n",
+    "    ylabel=\"$Z$ [m]\",\n",
+    "    grid=True,\n",
+    "    xlocator=1,\n",
+    "    ylocator=1,\n",
+    "    tickminor=True,\n",
+    ")"
+   ]
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "docs",
+   "language": "python",
+   "name": "python3"
+  },
+  "language_info": {
+   "codemirror_mode": {
+    "name": "ipython",
+    "version": 3
+   },
+   "file_extension": ".py",
+   "mimetype": "text/x-python",
+   "name": "python",
+   "nbconvert_exporter": "python",
+   "pygments_lexer": "ipython3",
+   "version": "3.14.6"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 5
+}

From 3c9d72940a6a16133256113ef8cbda257d2f3adb Mon Sep 17 00:00:00 2001
From: munechika-koyo <munechika.koyo@gmail.com>
Date: Tue, 11 Aug 2026 11:58:55 +0200
Subject: [PATCH 28/36] =?UTF-8?q?=F0=9F=94=A7=20Rearrange=20typos=20job=20?=
 =?UTF-8?q?in=20pre-commit=20hooks=20after=20nbstripout?=
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit

---
 .lefthook.yaml | 8 ++++----
 1 file changed, 4 insertions(+), 4 deletions(-)

diff --git a/.lefthook.yaml b/.lefthook.yaml
index 789d49d..696776a 100644
--- a/.lefthook.yaml
+++ b/.lefthook.yaml
@@ -47,10 +47,6 @@ pre-commit:
         - run: test ${CI:-false} = true
       run: pixi {run} pyrefly-check {staged_files}
 
-    - name: typos
-      stage_fixed: true
-      run: pixi {run} typos {staged_files}
-
     - name: actionlint
       glob: ".github/workflows/*.{yaml,yml}"
       run: pixi {run} actionlint {staged_files}
@@ -80,3 +76,7 @@ pre-commit:
       only:
         - run: test ${CI:-false} = true
       run: pixi {run} nbstripout --verify {staged_files}
+
+    - name: typos
+      stage_fixed: true
+      run: pixi {run} typos {staged_files}

From 5e86f1f5c18224bd14853cdfe6ffdc9a3875fd56 Mon Sep 17 00:00:00 2001
From: munechika-koyo <munechika.koyo@gmail.com>
Date: Tue, 11 Aug 2026 14:03:52 +0200
Subject: [PATCH 29/36] Configure Dependabot for GitHub Actions

---
 .github/dependabot.yml | 11 +++++++++++
 1 file changed, 11 insertions(+)
 create mode 100644 .github/dependabot.yml

diff --git a/.github/dependabot.yml b/.github/dependabot.yml
new file mode 100644
index 0000000..e38985c
--- /dev/null
+++ b/.github/dependabot.yml
@@ -0,0 +1,11 @@
+version: 2
+
+updates:
+  - package-ecosystem: github-actions
+    directory: /
+    schedule:
+      interval: weekly
+    groups:
+      github-actions:
+        patterns:
+          - "*"

From caac0059e8e76981e0d5ab796bebe7d1407820d1 Mon Sep 17 00:00:00 2001
From: munechika-koyo <munechika.koyo@gmail.com>
Date: Tue, 11 Aug 2026 14:11:26 +0200
Subject: [PATCH 30/36] Update GitHub Actions versions

---
 .github/workflows/ci.yaml         | 10 +++++-----
 .github/workflows/deploy-pypi.yml | 22 +++++++++++-----------
 .github/workflows/docs.yml        |  8 ++++----
 .github/workflows/release.yml     |  4 ++--
 4 files changed, 22 insertions(+), 22 deletions(-)

diff --git a/.github/workflows/ci.yaml b/.github/workflows/ci.yaml
index 6431e67..3a16b6a 100644
--- a/.github/workflows/ci.yaml
+++ b/.github/workflows/ci.yaml
@@ -12,10 +12,10 @@ jobs:
     runs-on: ubuntu-latest
     steps:
       - name: Check out repo
-        uses: actions/checkout@v5
+        uses: actions/checkout@v7
 
       - name: 🟨 Set up Pixi
-        uses: prefix-dev/setup-pixi@v0.9.4
+        uses: prefix-dev/setup-pixi@v0.10.1
         with:
           environments: lint
 
@@ -32,10 +32,10 @@ jobs:
 
     steps:
       - name: check out repo
-        uses: actions/checkout@v5
+        uses: actions/checkout@v7
 
       - name: 🟨 Set up Pixi
-        uses: prefix-dev/setup-pixi@v0.9.4
+        uses: prefix-dev/setup-pixi@v0.10.1
         with:
           environments: ${{ matrix.environment }}
 
@@ -48,6 +48,6 @@ jobs:
         run: mv .coverage ".coverage.${{ matrix.environment }}.${{ matrix.os }}.xml"
 
       - name: 📤 Upload coverage to Codecov
-        uses: codecov/codecov-action@v5
+        uses: codecov/codecov-action@v7
         with:
           token: ${{ secrets.CODECOV_TOKEN }}
diff --git a/.github/workflows/deploy-pypi.yml b/.github/workflows/deploy-pypi.yml
index cc901b6..e9b4969 100644
--- a/.github/workflows/deploy-pypi.yml
+++ b/.github/workflows/deploy-pypi.yml
@@ -29,7 +29,7 @@ jobs:
     steps:
       - name: Checkout LLVM on macOS
         if: runner.os == 'macOS'
-        uses: actions/checkout@v4
+        uses: actions/checkout@v7
         with:
           repository: llvm/llvm-project
           ref: release/19.x
@@ -56,14 +56,14 @@ jobs:
           cmake --install build
 
       - name: Checkout repository
-        uses: actions/checkout@v5
+        uses: actions/checkout@v7
 
       - name: Install uv
-        uses: astral-sh/setup-uv@v7
+        uses: astral-sh/setup-uv@v9
 
       # For aarch64 support
       # https://cibuildwheel.pypa.io/en/stable/faq/#emulation
-      - uses: docker/setup-qemu-action@v3
+      - uses: docker/setup-qemu-action@v4
         with:
           platforms: all
         if: runner.os == 'Linux' && matrix.arch == 'aarch64'
@@ -82,7 +82,7 @@ jobs:
 
       - name: Build wheels (at Pull Request)
         if: github.event_name == 'pull_request'
-        uses: pypa/cibuildwheel@v3.2.1
+        uses: pypa/cibuildwheel@v4.2
         env:
           CIBW_ARCHS: ${{ matrix.arch }}
           CIBW_ENVIRONMENT_MACOS: MACOSX_DEPLOYMENT_TARGET=${{ matrix.macos }}
@@ -90,13 +90,13 @@ jobs:
 
       - name: Build wheels
         if: github.event_name != 'pull_request'
-        uses: pypa/cibuildwheel@v3.2.1
+        uses: pypa/cibuildwheel@v4.2
         env:
           CIBW_ARCHS: ${{ matrix.arch }}
           CIBW_ENVIRONMENT_MACOS: MACOSX_DEPLOYMENT_TARGET=${{ matrix.macos }}
 
       - name: Upload wheels
-        uses: actions/upload-artifact@v5
+        uses: actions/upload-artifact@v7
         with:
           name: cibw-wheels-${{ matrix.os }}-${{ strategy.job-index }}
           path: ./wheelhouse/*.whl
@@ -105,12 +105,12 @@ jobs:
     name: Build source distribution
     runs-on: ubuntu-latest
     steps:
-      - uses: actions/checkout@v5
+      - uses: actions/checkout@v7
 
       - name: Build sdist
         run: pipx run build --sdist
 
-      - uses: actions/upload-artifact@v5
+      - uses: actions/upload-artifact@v7
         with:
           name: cibw-sdist
           path: dist/*.tar.gz
@@ -120,7 +120,7 @@ jobs:
     name: Show artifacts
     runs-on: ubuntu-latest
     steps:
-      - uses: actions/download-artifact@v6
+      - uses: actions/download-artifact@v8
         with:
           pattern: cibw-*
           path: ${{ github.workspace }}/dist
@@ -144,7 +144,7 @@ jobs:
     # if: github.event_name == 'push' && startsWith(github.ref, 'refs/tags/v')
     steps:
       - name: Download artifacts
-        uses: actions/download-artifact@v6
+        uses: actions/download-artifact@v8
         with:
           # unpacks all CIBW artifacts into dist/
           pattern: cibw-*
diff --git a/.github/workflows/docs.yml b/.github/workflows/docs.yml
index 1139185..2ff9a49 100644
--- a/.github/workflows/docs.yml
+++ b/.github/workflows/docs.yml
@@ -20,10 +20,10 @@ jobs:
     runs-on: ubuntu-latest
     steps:
       - name: Check out repo
-        uses: actions/checkout@v5
+        uses: actions/checkout@v7
 
       - name: 🟨 Set up Pixi
-        uses: prefix-dev/setup-pixi@v0.9.4
+        uses: prefix-dev/setup-pixi@v0.10.1
         with:
           environments: docs
 
@@ -31,7 +31,7 @@ jobs:
         run: pixi run -e docs doc-build
 
       - name: 📦 Upload artifact
-        uses: actions/upload-pages-artifact@v4
+        uses: actions/upload-pages-artifact@v5
         with:
           path: docs/build/html
 
@@ -45,4 +45,4 @@ jobs:
     steps:
       - name: 🚀 Deploy to GitHub Pages
         id: deployment
-        uses: actions/deploy-pages@v4
+        uses: actions/deploy-pages@v5
diff --git a/.github/workflows/release.yml b/.github/workflows/release.yml
index 3515032..c9b58e6 100644
--- a/.github/workflows/release.yml
+++ b/.github/workflows/release.yml
@@ -14,7 +14,7 @@ jobs:
     runs-on: ubuntu-latest
     steps:
       - name: Checkout Repository
-        uses: actions/checkout@v5
+        uses: actions/checkout@v7
 
       - name: Extract release notes
         id: notes
@@ -37,7 +37,7 @@ jobs:
           } >> "$GITHUB_OUTPUT"
 
       - name: Create GitHub release
-        uses: softprops/action-gh-release@v2
+        uses: softprops/action-gh-release@v3
         with:
           tag_name: ${{ github.ref_name }}
           name: ${{ github.ref_name }}

From 3161b0f5c1ed3525a81d4b53cd157077db7bb288 Mon Sep 17 00:00:00 2001
From: munechika-koyo <munechika.koyo@gmail.com>
Date: Tue, 11 Aug 2026 14:15:31 +0200
Subject: [PATCH 31/36] =?UTF-8?q?=F0=9F=94=A7=20Update=20setup-uv=20action?=
 =?UTF-8?q?=20to=20version=209.0.0?=
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit

---
 .github/workflows/deploy-pypi.yml | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/.github/workflows/deploy-pypi.yml b/.github/workflows/deploy-pypi.yml
index e9b4969..6b88841 100644
--- a/.github/workflows/deploy-pypi.yml
+++ b/.github/workflows/deploy-pypi.yml
@@ -59,7 +59,7 @@ jobs:
         uses: actions/checkout@v7
 
       - name: Install uv
-        uses: astral-sh/setup-uv@v9
+        uses: astral-sh/setup-uv@v9.0.0
 
       # For aarch64 support
       # https://cibuildwheel.pypa.io/en/stable/faq/#emulation

From 0f0803c93e5d8b07c609055e295a8cc76a7ea0bb Mon Sep 17 00:00:00 2001
From: munechika-koyo <munechika.koyo@gmail.com>
Date: Tue, 11 Aug 2026 15:29:39 +0200
Subject: [PATCH 32/36] =?UTF-8?q?=F0=9F=94=A7=20Update=20deploy-pypi=20wor?=
 =?UTF-8?q?kflow=20and=20pyproject.toml=20for=20macOS=20OpenMP=20support?=
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit

---
 .github/workflows/deploy-pypi.yml | 17 ++++++++++++++++-
 pyproject.toml                    |  6 ++++++
 2 files changed, 22 insertions(+), 1 deletion(-)

diff --git a/.github/workflows/deploy-pypi.yml b/.github/workflows/deploy-pypi.yml
index 6b88841..344ca61 100644
--- a/.github/workflows/deploy-pypi.yml
+++ b/.github/workflows/deploy-pypi.yml
@@ -32,7 +32,8 @@ jobs:
         uses: actions/checkout@v7
         with:
           repository: llvm/llvm-project
-          ref: release/19.x
+          # Keep this pinned to a stable release for reproducible wheels.
+          ref: llvmorg-22.1.8
           path: llvm-project
 
       - name: Build OpenMP on macOS
@@ -49,11 +50,14 @@ jobs:
             -DCMAKE_INSTALL_PREFIX="${PREFIX}" \
             -DCMAKE_C_COMPILER=clang \
             -DCMAKE_CXX_COMPILER=clang++ \
+            -DCMAKE_OSX_ARCHITECTURES="${{ matrix.arch }}" \
+            -DLIBOMP_ENABLE_SHARED=ON \
             -DLIBOMP_INSTALL_ALIASES=OFF \
             -S openmp \
             -B build
           cmake --build build --parallel
           cmake --install build
+          echo "REPAIR_LIBRARY_PATH=${PREFIX}/lib" >> "${GITHUB_ENV}"
 
       - name: Checkout repository
         uses: actions/checkout@v7
@@ -95,6 +99,17 @@ jobs:
           CIBW_ARCHS: ${{ matrix.arch }}
           CIBW_ENVIRONMENT_MACOS: MACOSX_DEPLOYMENT_TARGET=${{ matrix.macos }}
 
+      - name: Verify OpenMP is bundled on macOS
+        if: runner.os == 'macOS'
+        shell: bash
+        run: |
+          for wheel in wheelhouse/*.whl; do
+            if ! unzip -l "${wheel}" | grep -Eq 'libomp[^/]*\.dylib$'; then
+              echo "::error::libomp.dylib is not bundled in ${wheel}"
+              exit 1
+            fi
+          done
+
       - name: Upload wheels
         uses: actions/upload-artifact@v7
         with:
diff --git a/pyproject.toml b/pyproject.toml
index 2a595d9..03d628e 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -96,6 +96,12 @@ build-frontend = "build[uv]"
 test-extras = ["test"]
 test-command = "pytest {project}/tests"
 
+[tool.cibuildwheel.macos]
+repair-wheel-command = """\
+DYLD_LIBRARY_PATH=$REPAIR_LIBRARY_PATH delocate-wheel \
+--require-archs {delocate_archs} -w {dest_dir} -v {wheel}\
+"""
+
 # -----------------------
 # === Cov/Test config ===
 # -----------------------

From 30ee51e0f94369868eebd9f4e24aa788f1f75651 Mon Sep 17 00:00:00 2001
From: munechika-koyo <munechika.koyo@gmail.com>
Date: Tue, 11 Aug 2026 16:20:39 +0200
Subject: [PATCH 33/36] =?UTF-8?q?=F0=9F=94=A7=20Update=20CIBW=5FSKIP=20set?=
 =?UTF-8?q?tings=20for=20compatibility=20with=20newer=20CPython=20versions?=
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit

---
 .github/workflows/deploy-pypi.yml | 2 +-
 pyproject.toml                    | 2 +-
 2 files changed, 2 insertions(+), 2 deletions(-)

diff --git a/.github/workflows/deploy-pypi.yml b/.github/workflows/deploy-pypi.yml
index 344ca61..92cd4b7 100644
--- a/.github/workflows/deploy-pypi.yml
+++ b/.github/workflows/deploy-pypi.yml
@@ -80,7 +80,7 @@ jobs:
         # These needs to rotate every new Python release.
         if: github.event_name == 'pull_request'
         run: |
-          CIBW_SKIP="pp* *-musllinux* cp311-* cp313-* cp314-* cp314t-*"
+          CIBW_SKIP="pp* *-musllinux* cp311* cp313* cp314* cp315*"
           echo "CIBW_SKIP=${CIBW_SKIP}" >> "${GITHUB_ENV}"
           echo "Setting CIBW_SKIP=${CIBW_SKIP}"
 
diff --git a/pyproject.toml b/pyproject.toml
index 03d628e..d9c4e53 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -91,7 +91,7 @@ define_macros = [["NPY_NO_DEPRECATED_API", "NPY_1_7_API_VERSION"]]
 targets = ["**/*.pyx"]
 
 [tool.cibuildwheel]
-skip = "pp* *_ppc64le *_i686 *_s390x *-musllinux* cp313-* cp314-* cp314t-*" # TODO: add cp313/314 when cherab supports it
+skip = "pp* *_ppc64le *_i686 *_s390x *-musllinux* cp313* cp314* cp315*" # TODO: Suppport latest CPython after cherab supports it
 build-frontend = "build[uv]"
 test-extras = ["test"]
 test-command = "pytest {project}/tests"

From 939db380cefc46e482694cb878f84f4e03488336 Mon Sep 17 00:00:00 2001
From: munechika-koyo <munechika.koyo@gmail.com>
Date: Tue, 11 Aug 2026 17:21:15 +0200
Subject: [PATCH 34/36] =?UTF-8?q?=F0=9F=90=9B=20Add=20missing=20marks=20to?=
 =?UTF-8?q?=20parameterized=20test=20for=20memory=20URI?=
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit

---
 tests/plasma/test_utility.py | 3 ++-
 1 file changed, 2 insertions(+), 1 deletion(-)

diff --git a/tests/plasma/test_utility.py b/tests/plasma/test_utility.py
index d2f74cb..b346f6a 100644
--- a/tests/plasma/test_utility.py
+++ b/tests/plasma/test_utility.py
@@ -10,7 +10,7 @@
 @pytest.mark.parametrize(
     ("constructor", "entry_kwargs", "expected"),
     [
-        (
+        pytest.param(
             "uri",
             {
                 "uri_builder": lambda _tmp_path, _path_iter_jintrac: (
@@ -19,6 +19,7 @@
                 "mode": "w",
             },
             "imas:memory?path=cherab_test_memory_uri",
+            marks=pytest.mark.requires_imas_memory_backend,
         ),
         (
             "uri",

From fae8815366c4a3c6a80ad4a75910861ab4b17ef7 Mon Sep 17 00:00:00 2001
From: munechika-koyo <munechika.koyo@gmail.com>
Date: Wed, 12 Aug 2026 00:51:47 +0200
Subject: [PATCH 35/36] =?UTF-8?q?=F0=9F=90=9B=20Fix=20exterior=20mesh=20re?=
 =?UTF-8?q?ndering=20in=20grid=203D=20notebook?=
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit

---
 docs/notebooks/radiation/grid_3d.ipynb | 93 +++++++++++++++++---------
 1 file changed, 61 insertions(+), 32 deletions(-)

diff --git a/docs/notebooks/radiation/grid_3d.ipynb b/docs/notebooks/radiation/grid_3d.ipynb
index 330d71b..57f62d1 100644
--- a/docs/notebooks/radiation/grid_3d.ipynb
+++ b/docs/notebooks/radiation/grid_3d.ipynb
@@ -151,11 +151,7 @@
    "cell_type": "code",
    "execution_count": null,
    "id": "11",
-   "metadata": {
-    "tags": [
-     "nbsphinx-thumbnail"
-    ]
-   },
+   "metadata": {},
    "outputs": [],
    "source": [
     "# Extract cell center points for one toroidal slice\n",
@@ -225,8 +221,11 @@
    "id": "14",
    "metadata": {},
    "source": [
-    "Visualize the 3D grid in a quarter of the torus.\n",
-    "[plotly](https://plotly.com/python/) is used to visualize the tetrahedral mesh in the notebook."
+    "## Visualize the exterior mesh\n",
+    "\n",
+    "To keep the visualization readable and efficient to render, only the exterior faces of the quarter-torus grid are shown. Faces shared by adjacent cells and all internal edges are omitted, while the edges of the exterior cells are retained to show the surface mesh structure.\n",
+    "\n",
+    "This filtering affects only the visualization; the underlying volume mesh is not modified or downsampled. [Plotly](https://plotly.com/python/) is used to render the exterior mesh."
    ]
   },
   {
@@ -246,32 +245,62 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "tetra = grid_cut.tetrahedra\n",
-    "faces = np.vstack(\n",
+    "verts = grid_cut.vertices\n",
+    "cells = grid_cut.cells\n",
+    "\n",
+    "# Extract the six quadrilateral faces of each hexahedral cell.  Working with\n",
+    "# cell faces, rather than tetrahedral faces, allows shared internal faces to\n",
+    "# cancel even when adjacent cells use different tetrahedral diagonals.\n",
+    "cell_face_indices = np.array(\n",
     "    [\n",
-    "        tetra[:, [0, 1, 2]],\n",
-    "        tetra[:, [0, 1, 3]],\n",
-    "        tetra[:, [0, 2, 3]],\n",
-    "        tetra[:, [1, 2, 3]],\n",
+    "        [0, 1, 2, 3],\n",
+    "        [4, 5, 6, 7],\n",
+    "        [0, 1, 5, 4],\n",
+    "        [1, 2, 6, 5],\n",
+    "        [2, 3, 7, 6],\n",
+    "        [3, 0, 4, 7],\n",
     "    ]\n",
     ")\n",
+    "cell_faces = cells[:, cell_face_indices].reshape(-1, 4)\n",
     "\n",
-    "# Keep only boundary faces (faces appearing exactly once).\n",
-    "faces_sorted = np.sort(faces, axis=1)\n",
-    "_, inverse, counts = np.unique(faces_sorted, axis=0, return_inverse=True, return_counts=True)\n",
-    "surface_faces = faces[counts[inverse] == 1]\n",
+    "# A face that occurs only once belongs to the exterior of the selected grid.\n",
+    "_, face_indices, counts = np.unique(\n",
+    "    np.sort(cell_faces, axis=1),\n",
+    "    axis=0,\n",
+    "    return_index=True,\n",
+    "    return_counts=True,\n",
+    ")\n",
+    "surface_quads = cell_faces[face_indices[counts == 1]]\n",
     "\n",
-    "# Build unique boundary edges from boundary triangles.\n",
-    "tri_edges = np.vstack(\n",
+    "# Split exterior quads into triangles for Plotly's Mesh3d and discard\n",
+    "# zero-area triangles caused by padded triangular cells.\n",
+    "surface_faces = np.vstack([surface_quads[:, [0, 1, 2]], surface_quads[:, [0, 2, 3]]])\n",
+    "face_vertices = verts[surface_faces]\n",
+    "face_area_twice = np.linalg.norm(\n",
+    "    np.cross(\n",
+    "        face_vertices[:, 1] - face_vertices[:, 0],\n",
+    "        face_vertices[:, 2] - face_vertices[:, 0],\n",
+    "    ),\n",
+    "    axis=1,\n",
+    ")\n",
+    "surface_faces = surface_faces[face_area_twice > 1e-12]\n",
+    "\n",
+    "# Draw only exterior cell edges.  These reveal the original mesh structure\n",
+    "# without sending the hidden volume edges to the renderer.\n",
+    "surface_edges = np.vstack(\n",
     "    [\n",
-    "        surface_faces[:, [0, 1]],\n",
-    "        surface_faces[:, [1, 2]],\n",
-    "        surface_faces[:, [2, 0]],\n",
+    "        surface_quads[:, [0, 1]],\n",
+    "        surface_quads[:, [1, 2]],\n",
+    "        surface_quads[:, [2, 3]],\n",
+    "        surface_quads[:, [3, 0]],\n",
     "    ]\n",
     ")\n",
-    "boundary_edges = np.unique(np.sort(tri_edges, axis=1), axis=0)\n",
+    "edge_vertices = verts[surface_edges]\n",
+    "surface_edges = surface_edges[\n",
+    "    np.linalg.norm(edge_vertices[:, 1] - edge_vertices[:, 0], axis=1) > 1e-12\n",
+    "]\n",
+    "surface_edges = np.unique(np.sort(surface_edges, axis=1), axis=0)\n",
     "\n",
-    "verts = grid_cut.vertices\n",
     "x, y, z = verts[:, 0], verts[:, 1], verts[:, 2]\n",
     "\n",
     "# Customize edge appearance here.\n",
@@ -279,15 +308,15 @@
     "edge_width = 1.5\n",
     "\n",
     "# Convert edge index pairs to line segments separated by NaN for Plotly.\n",
-    "edge_xyz = verts[boundary_edges]\n",
+    "edge_xyz = verts[surface_edges]\n",
     "xe = np.column_stack(\n",
-    "    [edge_xyz[:, 0, 0], edge_xyz[:, 1, 0], np.full(len(boundary_edges), np.nan)]\n",
+    "    [edge_xyz[:, 0, 0], edge_xyz[:, 1, 0], np.full(len(surface_edges), np.nan)]\n",
     ").ravel()\n",
     "ye = np.column_stack(\n",
-    "    [edge_xyz[:, 0, 1], edge_xyz[:, 1, 1], np.full(len(boundary_edges), np.nan)]\n",
+    "    [edge_xyz[:, 0, 1], edge_xyz[:, 1, 1], np.full(len(surface_edges), np.nan)]\n",
     ").ravel()\n",
     "ze = np.column_stack(\n",
-    "    [edge_xyz[:, 0, 2], edge_xyz[:, 1, 2], np.full(len(boundary_edges), np.nan)]\n",
+    "    [edge_xyz[:, 0, 2], edge_xyz[:, 1, 2], np.full(len(surface_edges), np.nan)]\n",
     ").ravel()\n",
     "\n",
     "fig = go.Figure(\n",
@@ -311,16 +340,16 @@
     "            z=ze,\n",
     "            mode=\"lines\",\n",
     "            line=dict(color=edge_color, width=edge_width),\n",
-    "            name=\"Surface boundary\",\n",
+    "            name=\"Surface mesh\",\n",
     "            hoverinfo=\"skip\",\n",
     "        ),\n",
     "    ]\n",
     ")\n",
     "\n",
     "fig.update_layout(\n",
-    "    title=\"Quarter-torus tetrahedral mesh\",\n",
+    "    title=\"Quarter-torus surface mesh\",\n",
     "    width=500,\n",
-    "    height=500,\n",
+    "    height=600,\n",
     "    scene=dict(\n",
     "        xaxis_title=\"X [m]\",\n",
     "        yaxis_title=\"Y [m]\",\n",
@@ -329,7 +358,7 @@
     "    ),\n",
     "    margin=dict(l=0, r=0, b=0, t=30),\n",
     "    scene_camera=dict(\n",
-    "        eye=dict(x=0.0, y=-2.0, z=0.25),\n",
+    "        eye=dict(x=0.2, y=-2.25, z=0.5),\n",
     "    ),\n",
     "    template=\"plotly_dark\",\n",
     ")\n",

From 4696b420ecff3b273439f2e546cdee70a95e66de Mon Sep 17 00:00:00 2001
From: munechika-koyo <munechika.koyo@gmail.com>
Date: Wed, 12 Aug 2026 00:53:01 +0200
Subject: [PATCH 36/36] =?UTF-8?q?=F0=9F=94=A7=20Update=20changelog=20date?=
 =?UTF-8?q?=20for=20version=200.6.0?=
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit

---
 CHANGELOG.md | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/CHANGELOG.md b/CHANGELOG.md
index 2e79ac1..1d1491f 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -5,7 +5,7 @@ All notable changes to this project will be documented in this file.
 The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.1.0/),
 and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0.html).
 
-## [0.6.0]
+## [0.6.0] - 2026-08-12
 
 ### Added