r"""
Routines to find and classify critical points of 2D spline fields.
A critical point is a location where the gradient of the field is zero.
For a poloidal flux function $\psi(R, Z)$ these correspond to:
- O-points: local extrema of $\psi$, i.e. the magnetic axis.
- X-points: saddle points of $\psi$, i.e. the separatrix / divertor X-point.
Two approaches are provided:
``find_extremum``
Finds a local minimum or maximum of any 2D spline field by minimizing
the field value directly. The recommended approach for finding X-points
is to pass the poloidal field magnitude interpolator with ``mode="min"``,
since $|B_\mathrm{pol}| = 0$ exactly at the X-point making it a
well-defined minimum. For O-points, pass $\psi$ with ``mode="auto"``
to detect the sign convention automatically.
``find_critical_point``
Finds any critical point of $\psi$ (O- or X-point) by minimizing
$|\nabla\psi|^2$, which is zero at any critical point regardless of its
character. Use this as a fallback when a poloidal field interpolator is
not available.
Typical usage
-------------
Example usage::
# X-point via poloidal field (preferred)
x_point = find_extremum(bpol_interpolator, r0, z0, mode="min")
# X-point via psi directly (fallback)
x_point = find_critical_point(psi_interpolator, r0, z0)
# O-point -- auto detects sign convention from Laplacian
o_point = find_extremum(psi_interpolator, r0, z0, mode="auto")
# O-point -- explicit if sign convention is known
o_point = find_extremum(psi_interpolator, r0, z0, mode="max")
"""
from __future__ import annotations
from dataclasses import dataclass
from enum import Enum
import numpy as np
from scipy.interpolate import RectBivariateSpline
from scipy.optimize import minimize
[docs]
class CriticalPointType(Enum):
"""Enumeration of critical point types based on Hessian character."""
MINIMUM = "minimum"
MAXIMUM = "maximum"
SADDLE = "saddle"
DEGENERATE = "degenerate"
UNCERTAIN = "uncertain"
[docs]
@dataclass
class CriticalPoint:
"""A located and classified critical point of a 2D scalar field."""
r: float
z: float
point_type: CriticalPointType
def __repr__(self) -> str:
"""Return a string representation of the critical point."""
return (
f"CriticalPoint({self.point_type.value} "
f"at R={self.r:.6f}, Z={self.z:.6f})"
)
_DEGENERATE_THRESHOLD = 1e-10
# Fraction of Hessian scale below which X-point is uncertain
_UNCERTAIN_RTOL = 1e-4
_VALID_MODES = ("min", "max", "auto")
def _classify(spline: RectBivariateSpline, r: float, z: float) -> CriticalPointType:
r"""
Classify a critical point by the Hessian determinant and trace.
- $\det(H) > 0$, $\mathrm{tr}(H) > 0$: minimum.
- $\det(H) > 0$, $\mathrm{tr}(H) < 0$: maximum.
- $\det(H) < -\epsilon |\partial^2 f/\partial R^2 \cdot \partial^2 f/\partial Z^2|$: saddle.
- $|\det(H)| < $``_DEGENERATE_THRESHOLD``: degenerate.
- small negative $\det(H)$: uncertain.
"""
d2f_dR2 = spline(r, z, dx=2, dy=0, grid=False).item()
d2f_dZ2 = spline(r, z, dx=0, dy=2, grid=False).item()
d2f_dRdZ = spline(r, z, dx=1, dy=1, grid=False).item()
det_H = d2f_dR2 * d2f_dZ2 - d2f_dRdZ**2
trace_H = d2f_dR2 + d2f_dZ2
if abs(det_H) < _DEGENERATE_THRESHOLD:
return CriticalPointType.DEGENERATE
if det_H > 0:
return CriticalPointType.MINIMUM if trace_H > 0 else CriticalPointType.MAXIMUM
hessian_scale = abs(d2f_dR2 * d2f_dZ2)
if hessian_scale > 0 and abs(det_H) < _UNCERTAIN_RTOL * hessian_scale:
return CriticalPointType.UNCERTAIN
return CriticalPointType.SADDLE
def _detect_mode(field_interpolator: RectBivariateSpline, r: float, z: float) -> str:
r"""
Detect whether to minimize or maximize by checking the Laplacian sign.
Assumes the field is $\psi$ -- do not use with $|B_\mathrm{pol}|$ or
other fields.
- Laplacian > 0: local minimum -- returns ``"min"``.
- Laplacian < 0: local maximum -- returns ``"max"``.
"""
laplacian = (
field_interpolator(r, z, dx=2, dy=0, grid=False).item() +
field_interpolator(r, z, dx=0, dy=2, grid=False).item()
)
return "min" if laplacian > 0 else "max"
def _attempt_extremum(
field_interpolator: RectBivariateSpline,
initial_r: float,
initial_z: float,
mode: str,
tolerance: float,
**scipy_kwargs,
) -> CriticalPointLocation | None:
"""
Attempt to find an extremum in the given mode.
Returns a ``CriticalPointLocation`` or ``None`` on failure.
"""
sign = 1.0 if mode == "min" else -1.0
def objective(pt):
return sign * field_interpolator(pt[0], pt[1], grid=False).item()
def jacobian(pt):
df_dR = field_interpolator(pt[0], pt[1], dx=1, dy=0, grid=False).item()
df_dZ = field_interpolator(pt[0], pt[1], dx=0, dy=1, grid=False).item()
return sign * np.array([df_dR, df_dZ])
result = minimize(
objective,
x0=np.array([float(initial_r), float(initial_z)]),
jac=jacobian,
tol=tolerance,
**scipy_kwargs,
)
if not result.success:
return None
r, z = float(result.x[0]), float(result.x[1])
return CriticalPoint(r=r, z=z, point_type=_classify(field_interpolator, r, z))
[docs]
def find_extremum(
field_interpolator: RectBivariateSpline,
initial_r: float,
initial_z: float,
mode: str = "min",
tolerance: float = 1e-6,
**scipy_kwargs,
) -> CriticalPointLocation | None:
r"""
Find a local extremum of a 2D spline field.
Parameters
----------
field_interpolator:
A ``RectBivariateSpline`` representing any scalar field $f(R, Z)$.
initial_r:
Starting guess for the $R$ coordinate of the extremum.
initial_z:
Starting guess for the $Z$ coordinate of the extremum.
mode:
One of ``"min"``, ``"max"``, or ``"auto"``.
- ``"min"`` -- find the nearest local minimum.
- ``"max"`` -- find the nearest local maximum.
- ``"auto"`` -- infer min or max from the Laplacian sign at the
starting point, then try the other if the first attempt fails.
Only valid when ``field_interpolator`` is $\psi$.
tolerance:
Convergence tolerance passed to ``scipy.optimize.minimize``.
**scipy_kwargs:
Passed directly to ``scipy.optimize.minimize``.
Returns
-------
CriticalPointLocation or None
Located $(R, Z)$, or ``None`` on failure.
"""
assert mode in _VALID_MODES, (
f"mode must be one of {_VALID_MODES}, got '{mode}'"
)
if mode in ("min", "max"):
location = _attempt_extremum(
field_interpolator, initial_r, initial_z,
mode=mode, tolerance=tolerance, **scipy_kwargs,
)
if location is None:
print(
f"find_extremum: optimization failed starting from "
f"(R={initial_r:.4f}, Z={initial_z:.4f}) with mode='{mode}'."
)
return location
primary = _detect_mode(field_interpolator, initial_r, initial_z)
fallback = "max" if primary == "min" else "min"
location = _attempt_extremum(
field_interpolator, initial_r, initial_z,
mode=primary, tolerance=tolerance, **scipy_kwargs,
)
if location is not None:
return location
print(
f"find_extremum: auto mode with mode='{primary}' failed "
f"starting from (R={initial_r:.4f}, Z={initial_z:.4f}) "
f"-- retrying with mode='{fallback}'."
)
location = _attempt_extremum(
field_interpolator, initial_r, initial_z,
mode=fallback, tolerance=tolerance, **scipy_kwargs,
)
if location is not None:
return location
print(
f"find_extremum: auto mode failed both '{primary}' and '{fallback}' "
f"starting from (R={initial_r:.4f}, Z={initial_z:.4f})."
)
return None
[docs]
def find_critical_point(
psi_interpolator: RectBivariateSpline,
initial_r: float,
initial_z: float,
tolerance: float = 1e-6,
**scipy_kwargs,
) -> CriticalPoint | None:
r"""
Find a critical point of $\psi(R, Z)$ and classify it.
Minimizes $|\nabla\psi|^2 = (\partial\psi/\partial R)^2 +
(\partial\psi/\partial Z)^2$, which is zero at any critical point of
$\psi$ regardless of its character (minimum, maximum, or saddle). The
character is then determined from the sign of the Hessian determinant.
Use this as a fallback when a poloidal field magnitude interpolator is
not available. When $|B_\mathrm{pol}|$ is available, prefer
``find_extremum`` with ``mode="min"``, which is more numerically robust
for X-point finding because the objective is physically motivated and
strictly non-negative with a clear zero.
Parameters
----------
psi_interpolator:
A ``RectBivariateSpline`` representing the poloidal flux $\psi(R, Z)$.
initial_r:
Starting guess for the $R$ coordinate of the critical point.
initial_z:
Starting guess for the $Z$ coordinate of the critical point.
tolerance:
Convergence tolerance passed to ``scipy.optimize.minimize``.
**scipy_kwargs:
Passed directly to ``scipy.optimize.minimize``.
Returns
-------
CriticalPoint or None
Located $(R, Z)$ and classification, or ``None`` on failure.
"""
def objective(pt):
dpsi_dR = psi_interpolator(pt[0], pt[1], dx=1, dy=0, grid=False).item()
dpsi_dZ = psi_interpolator(pt[0], pt[1], dx=0, dy=1, grid=False).item()
return dpsi_dR**2 + dpsi_dZ**2
def jacobian(pt):
dpsi_dR = psi_interpolator(pt[0], pt[1], dx=1, dy=0, grid=False).item()
dpsi_dZ = psi_interpolator(pt[0], pt[1], dx=0, dy=1, grid=False).item()
d2psi_dR2 = psi_interpolator(pt[0], pt[1], dx=2, dy=0, grid=False).item()
d2psi_dZ2 = psi_interpolator(pt[0], pt[1], dx=0, dy=2, grid=False).item()
d2psi_dRdZ = psi_interpolator(pt[0], pt[1], dx=1, dy=1, grid=False).item()
df_dR = 2 * (dpsi_dR * d2psi_dR2 + dpsi_dZ * d2psi_dRdZ)
df_dZ = 2 * (dpsi_dR * d2psi_dRdZ + dpsi_dZ * d2psi_dZ2)
return np.array([df_dR, df_dZ])
result = minimize(
objective,
x0=np.array([float(initial_r), float(initial_z)]),
jac=jacobian,
tol=tolerance,
**scipy_kwargs,
)
if not result.success:
print(
f"find_critical_point: optimization failed starting from "
f"(R={initial_r:.4f}, Z={initial_z:.4f}).\n"
f" Reason : {result.message}\n"
f" |grad psi|^2 : {result.fun:.3e}"
)
return None
r, z = float(result.x[0]), float(result.x[1])
point_type = _classify(psi_interpolator, r, z)
if point_type == CriticalPointType.DEGENERATE:
print(
f"find_critical_point: point at (R={r:.4f}, Z={z:.4f}) "
f"has |det(H)| < {_DEGENERATE_THRESHOLD:.1e} -- "
f"classification is unreliable."
)
return CriticalPoint(r=r, z=z, point_type=point_type)
