fusion_core/
xpoint.rs

1// SPDX-License-Identifier: AGPL-3.0-or-later | Commercial license available
2// © Concepts 1996–2026 Miroslav Šotek. All rights reserved.
3// © Code 2020–2026 Miroslav Šotek. All rights reserved.
4// ORCID: 0009-0009-3560-0851
5// Contact: www.anulum.li | protoscience@anulum.li
6// SCPN Fusion Core — Xpoint
7//! X-point (magnetic null) detection.
8//!
9//! Port of fusion_kernel.py `find_x_point()` (lines 95-116).
10//! Locates the saddle point where B=0 in the divertor region.
11
12use fusion_math::interp::gradient_2d;
13use fusion_types::error::{FusionError, FusionResult};
14use fusion_types::state::Grid2D;
15use ndarray::Array2;
16
17/// Find the X-point (magnetic null) in the lower divertor region.
18///
19/// Algorithm:
20/// 1. Compute gradient magnitude |∇Ψ| using central differences
21/// 2. Mask upper half (only search Z < Z_min * 0.5)
22/// 3. Find argmin of |∇Ψ| in masked region
23///
24/// Returns `((R, Z), psi_value)` at the X-point.
25///
26/// NOTE: The Python code calls `np.gradient(Psi, self.dR, self.dZ)` which returns
27/// gradients along axis 0 (with spacing dR) and axis 1 (with spacing dZ).
28/// The Python variable names are swapped (dPsi_dR is actually along Z-axis,
29/// dPsi_dZ along R-axis), but the gradient magnitude is the same either way.
30/// In Rust we use correct naming via gradient_2d().
31pub fn find_x_point(
32    psi: &Array2<f64>,
33    grid: &Grid2D,
34    z_min: f64,
35) -> FusionResult<((f64, f64), f64)> {
36    if !z_min.is_finite() {
37        return Err(FusionError::ConfigError(format!(
38            "x-point z_min must be finite, got {z_min}"
39        )));
40    }
41    if grid.nz < 2 || grid.nr < 2 {
42        return Err(FusionError::ConfigError(format!(
43            "x-point grid requires nz,nr >= 2, got nz={} nr={}",
44            grid.nz, grid.nr
45        )));
46    }
47    if psi.nrows() != grid.nz || psi.ncols() != grid.nr {
48        return Err(FusionError::ConfigError(format!(
49            "x-point psi shape mismatch: expected ({}, {}), got ({}, {})",
50            grid.nz,
51            grid.nr,
52            psi.nrows(),
53            psi.ncols()
54        )));
55    }
56    if psi.iter().any(|v| !v.is_finite()) {
57        return Err(FusionError::ConfigError(
58            "x-point psi contains non-finite values".to_string(),
59        ));
60    }
61    if grid.z.iter().any(|v| !v.is_finite()) || grid.r.iter().any(|v| !v.is_finite()) {
62        return Err(FusionError::ConfigError(
63            "x-point grid axes must be finite".to_string(),
64        ));
65    }
66
67    // Compute gradient components (correct naming)
68    let (dpsi_dz, dpsi_dr) = gradient_2d(psi, grid);
69
70    let nz = grid.nz;
71    let nr = grid.nr;
72
73    // Gradient magnitude |∇Ψ|
74    // Note: In Python the variable names are swapped but |∇Ψ| = sqrt(a² + b²) is the same
75    let mut min_b_mag = f64::MAX;
76    let mut best_iz = 0;
77    let mut best_ir = 0;
78    let mut found = false;
79
80    let z_threshold = z_min * 0.5;
81
82    for iz in 0..nz {
83        for ir in 0..nr {
84            let z = grid.zz[[iz, ir]];
85
86            // Only search in divertor region (Z < Z_min * 0.5)
87            if z < z_threshold {
88                let b_mag_sq = dpsi_dz[[iz, ir]].powi(2) + dpsi_dr[[iz, ir]].powi(2);
89                if !b_mag_sq.is_finite() || b_mag_sq < 0.0 {
90                    return Err(FusionError::ConfigError(format!(
91                        "x-point gradient magnitude became invalid at ({iz}, {ir})"
92                    )));
93                }
94                let b_mag = b_mag_sq.sqrt();
95
96                if b_mag < min_b_mag {
97                    min_b_mag = b_mag;
98                    best_iz = iz;
99                    best_ir = ir;
100                    found = true;
101                }
102            }
103        }
104    }
105
106    if !found {
107        return Err(FusionError::ConfigError(
108            "x-point search region is empty for provided z_min threshold".to_string(),
109        ));
110    }
111
112    let r = grid.r[best_ir];
113    let z = grid.z[best_iz];
114    let psi_x = psi[[best_iz, best_ir]];
115    if !r.is_finite() || !z.is_finite() || !psi_x.is_finite() {
116        return Err(FusionError::ConfigError(
117            "x-point output contains non-finite values".to_string(),
118        ));
119    }
120    Ok(((r, z), psi_x))
121}
122
123#[cfg(test)]
124mod tests {
125    use super::*;
126
127    #[test]
128    fn test_find_x_point_basic() {
129        let grid = Grid2D::new(33, 33, 1.0, 9.0, -5.0, 5.0);
130        // Create a simple field with a null at (5.0, -2.5)
131        let psi = Array2::from_shape_fn((33, 33), |(iz, ir)| {
132            let r = grid.rr[[iz, ir]];
133            let z = grid.zz[[iz, ir]];
134            // Saddle point-like field
135            (r - 5.0).powi(2) - (z + 2.5).powi(2)
136        });
137
138        let ((r_x, z_x), _psi_x) = find_x_point(&psi, &grid, -5.0).expect("valid x-point inputs");
139
140        // X-point should be near (5.0, -2.5) — the saddle point
141        assert!(
142            (r_x - 5.0).abs() < 1.0,
143            "X-point R={r_x}, expected near 5.0"
144        );
145        assert!(
146            (z_x + 2.5).abs() < 1.0,
147            "X-point Z={z_x}, expected near -2.5"
148        );
149    }
150
151    #[test]
152    fn test_find_x_point_rejects_invalid_runtime_inputs() {
153        let grid = Grid2D::new(33, 33, 1.0, 9.0, -5.0, 5.0);
154        let psi = Array2::from_shape_fn((33, 33), |(iz, ir)| {
155            let r = grid.rr[[iz, ir]];
156            let z = grid.zz[[iz, ir]];
157            (r - 5.0).powi(2) - (z + 2.5).powi(2)
158        });
159
160        let err = find_x_point(&psi, &grid, f64::NAN).expect_err("non-finite z_min must fail");
161        assert!(matches!(err, FusionError::ConfigError(_)));
162
163        let err =
164            find_x_point(&psi, &grid, -1000.0).expect_err("empty divertor search region must fail");
165        assert!(matches!(err, FusionError::ConfigError(_)));
166    }
167}
fusion_core/xpoint.rs

fusion_core/
xpoint.rs
