Skip to content

feat(equilibrium): enhance contour_tree for advanced plasma topologies#243

Open
Simon-McIntosh wants to merge 4 commits into
iterorganization:developfrom
Simon-McIntosh:feature/contour-tree-enhancements
Open

feat(equilibrium): enhance contour_tree for advanced plasma topologies#243
Simon-McIntosh wants to merge 4 commits into
iterorganization:developfrom
Simon-McIntosh:feature/contour-tree-enhancements

Conversation

@Simon-McIntosh
Copy link
Copy Markdown
Contributor

@Simon-McIntosh Simon-McIntosh commented May 7, 2026

Summary

Enhances the equilibrium/time_slice/contour_tree schema to support double-null divertor (DND), doublet, and snowflake plasma topologies. The existing contour tree fields (node, edges, levelset) had minimal documentation and no mechanism to represent multi-region plasmas. This PR adds two new optional fields and substantially hardens the documentation of all existing fields.

Changes

New fields

Field Type Description
contour_tree/n_plasma_regions INT_0D Number of distinct confined plasma regions (1 for standard tokamak, >1 for doublets)
contour_tree/node/region INT_1D Region index(es) each critical point belongs to or borders; [-1] for vacuum nulls

Both fields are optional (minOccurs="0") — existing producers that only write standard single-region equilibria are unaffected.

Documentation hardening

node ordering — Tightened from advisory ("should") to mandatory ("must"):

  • Node 0: primary O-point (magnetic axis, minimum ψ)
  • Node 1: primary X-point (lowest ψ among X-points)
  • Nodes 2+: ascending ψ order
  • DND: secondary X-point must have higher node index than primary
  • Doublet: each confined region owns an O-point; lowest ψ O-point is primary

levelset segments — Replaced the single-sentence original with precise extraction rules:

  • Only topologically connected segments; disconnected contours at the same ψ must be excluded
  • Segment 0 is always the LCFS; legs ordered counterclockwise from outboard midplane
  • Defined θ=0 reference (R_max at X-point Z) and direction (counterclockwise in R,Z)
  • SND, DND, and snowflake segment counts documented explicitly
  • O-points: levelset must be empty (no contour passes through a flux extremum)

edges — Clarified as encoding the Reeb graph (contour tree) of the poloidal flux function. Added the mathematical definition: an edge (i, j) means sweeping ψ from node i to node j encounters no other critical point. Documented the forest structure for multi-region configurations. Edge index ordering upgraded from "should" to "must".

Disambiguation

  • Replaced ambiguous "nodes" (could mean grid points or contour tree nodes) with "critical points" vs "computational grid nodes" where context requires both
  • Inlined critical_type enum values (0=minimum, 1=saddle, 2=maximum) directly into node documentation for readability

Backward compatibility

  • No breaking changes. Both new fields are optional with minOccurs="0".
  • Existing producers continue to work — single-region equilibria need not populate region or n_plasma_regions.
  • Documentation changes only tighten language from advisory ("should") to mandatory ("must") for ordering conventions that all known implementations already follow.
  • The psi field documentation and its COCOS transformation annotations are unchanged.

Motivation

The EFIT (IO) equilibrium reconstruction code is implementing robust X-point detection and contour tree computation (NOVA-sourced algebraic topology classifier + Fortran merge-tree). The existing DD schema lacked the expressiveness to represent:

  1. DND configurations where two X-points define the separatrix
  2. TCV-style doublet plasmas with multiple confined regions sharing an X-point
  3. Snowflake configurations with higher-order null points

Without deterministic ordering rules and multi-region support, interoperability between codes producing and consuming contour trees is ambiguous.


📚 Documentation preview 📚: https://imas-data-dictionary--243.org.readthedocs.build/en/243/

Simon-McIntosh and others added 2 commits May 7, 2026 11:51
…pologies

Co-authored-by: Copilot <223556219+Copilot@users.noreply.github.com>
Tighten region, levelset, node, and edges documentation:
- region: clarify vacuum sentinel must be [-1], simplify region
  index assignment rule, fix single-region description
- levelset: remove redundant 'counterclockwise', remove
  'measure-zero' DND tangent, simplify departure angle spec
- node: disambiguate 'nodes' vs computational grid nodes,
  inline critical_type enum values, simplify no-X-point ordering
- edges: 'should' -> 'must' for index ordering convention,
  tighten forest description
@github-actions
Copy link
Copy Markdown

github-actions Bot commented May 7, 2026

Comment thread schemas/equilibrium/dd_equilibrium.xsd Outdated
<xs:element name="node" type="equilibrium_contour_tree_node" maxOccurs="unbounded">
<xs:annotation>
<xs:documentation>Set of nodes. A contour tree node is defined by its critical type and position within the poloidal plane. A critical type of 1 references an x-point whilst a critical type of 0 or 2 references an o-point. Both X-points and O-points are rarely coincident with nodes defining the poloidal upon which the poloidal flux map is defined. The order in which the critical points are stored in the nodes array of structure is only important for the primary plasma O-point and X-points. If present, the primary plasma O-point should be placed as the first element of the nodes array of structure. If present, the primary plasma X-point should be the second element in the nodes structure</xs:documentation>
<xs:documentation>Set of critical points of the poloidal flux function. Each node is defined by its critical type (0=minimum, 1=saddle/X-point, 2=maximum) and position within the poloidal plane. Critical points rarely coincide with the computational grid nodes on which the poloidal flux map is defined. Node 0 must contain the primary O-point (magnetic axis), which has the minimum psi value. Node 1 must contain the primary X-point, defined as the X-point with the lowest psi value (closest to the magnetic axis in psi). For double-null divertor (DND) configurations, the secondary X-point has a higher psi value and must be stored at a higher node index than the primary; the two psi values are always distinct in practice. For limiter-bounded plasmas where no X-point defines the LCFS, node 1 should contain the X-point nearest to the plasma boundary in psi, if any X-point exists. If no X-points exist, nodes are ordered by ascending psi value starting from node 1. Nodes 2 and above must contain additional critical points ordered by ascending psi value, i.e. from core outward. For doublet or multi-region plasmas, each confined plasma region has its own O-point; the O-point with the lowest psi value is primary (node 0).</xs:documentation>
Copy link
Copy Markdown
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Whether the value of poloidal flux at the magnetic axis (node 0) is a minimum or maximum will depend upon the direction of the plasma current. If you want to be specific on the type of extremum (min/max) then I think you will have to refer to normalised poloidal flux, psi_norm = (psi - psi_axis)/(psi_boundary - psi_axis).

Copy link
Copy Markdown
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

good catch, we have this allowance in the critical point type but the text needs to refer to normalized flux throughout. I will update.

@github-actions
Copy link
Copy Markdown

github-actions Bot commented May 7, 2026

@Simon-McIntosh Simon-McIntosh force-pushed the feature/contour-tree-enhancements branch from 3349651 to 07ab20c Compare May 7, 2026 12:28
@github-actions
Copy link
Copy Markdown

github-actions Bot commented May 7, 2026

…ntions

All references to psi magnitude ordering (lower/higher, ascending,
minimum/maximum) in the contour tree documentation now refer to
normalised poloidal flux psi_norm = (psi - psi_axis) / (psi_boundary
- psi_axis), which is always 0 at the magnetic axis and 1 at the
plasma boundary regardless of the sign of the plasma current.

Whether poloidal flux increases or decreases from the axis outward
depends on the sign of the plasma current, so poloidal flux
comparisons were ambiguous for specifying node ordering.

Changes:
- critical_type: clarify it classifies the poloidal flux; note
  that whether the axis is a minimum or maximum depends on Ip
- psi: drop Ip-specific directional comment; introduce psi_norm
  definition and state that all ordering rules use psi_norm
- node: replace all 'ascending psi' / 'lowest psi' ordering rules
  with psi_norm equivalents
- levelset: replace 'lower psi' with 'lower psi_norm' in DND text
- edges: clarify 'sweeping' language to be sign-independent
- contour_tree: fix typos (ether->either, extremum->extrema); add
  psi_norm definition and state ordering independence from Ip sign
@Simon-McIntosh Simon-McIntosh force-pushed the feature/contour-tree-enhancements branch from 07ab20c to 7176a58 Compare May 7, 2026 12:32
@github-actions
Copy link
Copy Markdown

github-actions Bot commented May 7, 2026

Comment thread schemas/equilibrium/dd_equilibrium.xsd Outdated
<xs:group ref="INT_1D"/>
</xs:complexType>
</xs:element>
<xs:element name="levelset" type="rz1d_dynamic_aos">
Copy link
Copy Markdown
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

If you want "levelset" to become an array of structure, (which I understand from the new definition proposed), it should be:

Suggested change
<xs:element name="levelset" type="rz1d_dynamic_aos">
<xs:element name="levelset" type="rz1d_dynamic_aos" maxOccurs="unbounded">

Some NBC metadata could also be inserted in the xs:appinfo section (although "contour_tree" is an alpha structure):
<change_nbc_version>4.2.0</change_nbc_version>
<change_nbc_description>type_changed</change_nbc_description>
<change_nbc_previous_type>rz1d_dynamic_aos</change_nbc_previous_type>

Copy link
Copy Markdown
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Good catch, thank you. Will add maxOccurs="unbounded" in the next push.

Copy link
Copy Markdown
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

N.B. this is a breaking change for codes, including the DD3to4 conversion in IMAS-Python and METIS (@jfartaud FYI)

Copy link
Copy Markdown
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I thought that we allowed these changes in alpha fields. Should we add the NBC metadata?

Copy link
Copy Markdown
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

It is allowed indeed, but please be aware that this will impact production code, which needs to be updated to support this change. Adding NBC metadata doesn't change that fact: the NBC metadata won't rewrite the source code of those programs 😉.

In the IMAS-Python case, we'd need to write something like:

if data_dictionary_version < Version("4.2.0"):
    levelset = node[1].levelset
else:  # levelset became an AoS in DD 4.2.0
    node[1].levelset.resize(1)
    levelset = node[1].levelset[0]
levelset.r = ...
levelset.z = ...

Comment thread schemas/equilibrium/dd_equilibrium.xsd Outdated
</xs:appinfo>
</xs:annotation>
<xs:sequence>
<xs:element name="n_plasma_regions" minOccurs="0">
Copy link
Copy Markdown
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The DD convention for numbers is to use "_n" as a suffix :

Suggested change
<xs:element name="n_plasma_regions" minOccurs="0">
<xs:element name="plasma_regions_n" minOccurs="0">

Copy link
Copy Markdown
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Thanks for the convention note. We are reworking how regions are handled as part of broader changes to support wall contact topology. This field may be removed or restructured.

Comment thread schemas/equilibrium/dd_equilibrium.xsd Outdated
<xs:sequence>
<xs:element name="n_plasma_regions" minOccurs="0">
<xs:annotation>
<xs:documentation>Number of distinct confined plasma regions in the equilibrium. For standard single-region tokamak configurations (limiter, SND, DND), n_plasma_regions=1. For doublet configurations (e.g. TCV doublet), n_plasma_regions equals the number of separate confined plasma domains, each containing at least one O-point. This value is consistent with the region indices assigned to the contour tree nodes.</xs:documentation>
Copy link
Copy Markdown
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Suggested change
<xs:documentation>Number of distinct confined plasma regions in the equilibrium. For standard single-region tokamak configurations (limiter, SND, DND), n_plasma_regions=1. For doublet configurations (e.g. TCV doublet), n_plasma_regions equals the number of separate confined plasma domains, each containing at least one O-point. This value is consistent with the region indices assigned to the contour tree nodes.</xs:documentation>
<xs:documentation>Number of distinct confined plasma regions in the equilibrium. For standard single-region tokamak configurations (limiter, SND, DND), plasma_regions_n=1. For doublet configurations (e.g. TCV doublet), plasma_regions_n equals the number of separate confined plasma domains, each containing at least one O-point. This value is consistent with the region indices assigned to the contour tree nodes.</xs:documentation>

Copy link
Copy Markdown
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

See above — this field is being reworked.

Comment thread schemas/equilibrium/dd_equilibrium.xsd Outdated
</xs:appinfo>
</xs:annotation>
<xs:sequence>
<xs:element name="n_plasma_regions" minOccurs="0">
Copy link
Copy Markdown
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I wonder if describing only the number of plasma regions is sufficient. Maybe it would be useful to have an array of structure identifier to describe the type of each region ?

Copy link
Copy Markdown
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

We are reworking region handling but a full identifier structure for region types is beyond the scope of this PR. The tree topology already captures the essential basin structure.

Copy link
Copy Markdown
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

We have adopted the identifier approach. A node_type identifier has been added to each node with values plasma (index 0) and vacuum (index 1), distinguishing confined plasma O-points from vacuum field extrema near external coils. The n_plasma_regions and region fields have been dropped in favour of this.

@github-actions
Copy link
Copy Markdown

@Simon-McIntosh Simon-McIntosh force-pushed the feature/contour-tree-enhancements branch from f1825b4 to 67e0421 Compare May 12, 2026 07:48
@github-actions
Copy link
Copy Markdown

…n fields

Add critical_type = 3 (limiter point) for points on the first wall
where the plasma boundary makes contact. Limiter point nodes are
inserted by splitting existing tree edges; they have degree 2 and
store the last closed flux surface in their levelset.

Limiter point validity is defined by the saddle-free subgraph: a
limiter point is valid only if it shares a connected component with
at least one O-point after removing all saddle nodes.

Add node_type identifier (plasma / vacuum) to distinguish confined
plasma O-points from vacuum field extrema near external coils.

Drop n_plasma_regions and region fields. The number of confined
plasma regions is the count of O-points with node_type = plasma.

Additional changes:
- Add maxOccurs=unbounded and coordinate1 to levelset
- Partitioned node ordering: critical points first, limiter points after
- Update psi, levelset, edges, and top-level documentation
- Add R, Z tie-break for equal psi_norm ordering
@Simon-McIntosh Simon-McIntosh force-pushed the feature/contour-tree-enhancements branch from 67e0421 to fff82aa Compare May 12, 2026 08:03
@github-actions
Copy link
Copy Markdown

Copy link
Copy Markdown
Collaborator

@SimonPinches SimonPinches left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

lgtm

<header>Node type for contour tree nodes, distinguishing confined plasma features from vacuum field features</header>
<dd_instance xpath="equilibrium/time_slice/contour_tree/node/node_type"/>

<int name="plasma" description="Magnetic axis of a confined plasma region" unique="yes">0</int>
Copy link
Copy Markdown
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The unique identifier metadata is deprecated, and should be removed

<dd_instance xpath="equilibrium/time_slice/contour_tree/node/node_type"/>

<int name="plasma" description="Magnetic axis of a confined plasma region" unique="yes">0</int>
<int name="vacuum" description="Vacuum field extremum, e.g. near external coils" unique="yes">1</int>
Copy link
Copy Markdown
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

So the primary X-point is a vacuum field feature ? I find it a bit counter-intuitive - we often speak of plasma X-point, and there is the SOL around it.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

None yet

Projects

None yet

Development

Successfully merging this pull request may close these issues.

7 participants