fusion_core/
xpoint.rs

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

fusion_core/
xpoint.rs
