fusion_core/
vacuum.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 — Vacuum
7//! Vacuum magnetic field from coils via toroidal Green's function.
8//!
9//! Port of fusion_kernel.py `calculate_vacuum_field()` (lines 59-93).
10//! Uses complete elliptic integrals K(m), E(m) for the exact toroidal
11//! single-coil flux function (Smythe/Jackson form).
12
13use fusion_math::elliptic::{ellipe, ellipk};
14use fusion_types::config::CoilConfig;
15use fusion_types::error::{FusionError, FusionResult};
16use fusion_types::state::Grid2D;
17use ndarray::Array2;
18
19/// Calculate vacuum magnetic flux from coils using the toroidal Green's function.
20///
21/// For each coil at `(Rc, Zc)` with current `I`, the flux contribution is:
22///
23///   Ψ = (μ₀ I / 2π) · √((R + Rc)² + (Z - Zc)²) · ((2 - k²)K(k²) - 2E(k²)) / k²
24///
25/// where k² = 4 R Rc / ((R + Rc)² + (Z - Zc)²), clipped to [1e-9, 0.999999].
26///
27/// Returns Ψ_vacuum `[nz, nr]`.
28pub fn calculate_vacuum_field(
29    grid: &Grid2D,
30    coils: &[CoilConfig],
31    mu0: f64,
32) -> FusionResult<Array2<f64>> {
33    if !mu0.is_finite() || mu0 <= 0.0 {
34        return Err(FusionError::ConfigError(format!(
35            "vacuum permeability mu0 must be finite and > 0, got {mu0}"
36        )));
37    }
38    if grid.nz == 0 || grid.nr == 0 {
39        return Err(FusionError::ConfigError(
40            "vacuum field grid must have nz,nr >= 1".to_string(),
41        ));
42    }
43    if grid.rr.nrows() != grid.nz || grid.rr.ncols() != grid.nr {
44        return Err(FusionError::ConfigError(format!(
45            "grid.rr shape mismatch: expected ({}, {}), got ({}, {})",
46            grid.nz,
47            grid.nr,
48            grid.rr.nrows(),
49            grid.rr.ncols()
50        )));
51    }
52    if grid.zz.nrows() != grid.nz || grid.zz.ncols() != grid.nr {
53        return Err(FusionError::ConfigError(format!(
54            "grid.zz shape mismatch: expected ({}, {}), got ({}, {})",
55            grid.nz,
56            grid.nr,
57            grid.zz.nrows(),
58            grid.zz.ncols()
59        )));
60    }
61    if grid.rr.iter().any(|v| !v.is_finite()) || grid.zz.iter().any(|v| !v.is_finite()) {
62        return Err(FusionError::ConfigError(
63            "vacuum field grid coordinates must be finite".to_string(),
64        ));
65    }
66
67    let nz = grid.nz;
68    let nr = grid.nr;
69    let mut psi_vac: Array2<f64> = Array2::zeros((nz, nr));
70
71    for coil in coils {
72        let rc = coil.r;
73        let zc = coil.z;
74        let current = coil.current;
75        if !rc.is_finite() || rc <= 0.0 {
76            return Err(FusionError::ConfigError(format!(
77                "coil radius must be finite and > 0, got {}",
78                coil.r
79            )));
80        }
81        if !zc.is_finite() || !current.is_finite() {
82            return Err(FusionError::ConfigError(
83                "coil z/current must be finite".to_string(),
84            ));
85        }
86
87        let prefactor = (mu0 * current) / (2.0 * std::f64::consts::PI);
88        if !prefactor.is_finite() {
89            return Err(FusionError::ConfigError(
90                "coil prefactor became non-finite".to_string(),
91            ));
92        }
93
94        for iz in 0..nz {
95            for ir in 0..nr {
96                let r = grid.rr[[iz, ir]];
97                let z = grid.zz[[iz, ir]];
98
99                let dz = z - zc;
100                let r_plus_rc = r + rc;
101                let r_plus_rc_sq = r_plus_rc * r_plus_rc;
102                let denom = r_plus_rc_sq + dz * dz;
103                if !denom.is_finite() || denom <= 0.0 {
104                    return Err(FusionError::ConfigError(
105                        "vacuum field denominator became non-finite or non-positive".to_string(),
106                    ));
107                }
108
109                // k² parameter
110                let k2_raw = (4.0 * r * rc) / denom;
111                if !k2_raw.is_finite() {
112                    return Err(FusionError::ConfigError(
113                        "vacuum field k^2 became non-finite".to_string(),
114                    ));
115                }
116                let mut k2 = k2_raw;
117                // CRITICAL: clip to avoid singularity at ellipk(1.0) and division by k2=0
118                k2 = k2.clamp(1e-9, 0.999999);
119
120                let k_val = ellipk(k2);
121                let e_val = ellipe(k2);
122                if !k_val.is_finite() || !e_val.is_finite() {
123                    return Err(FusionError::ConfigError(
124                        "vacuum field elliptic integrals became non-finite".to_string(),
125                    ));
126                }
127
128                let sqrt_term = (r * rc).sqrt();
129                let k = k2.sqrt();
130                let term = ((2.0 - k2) * k_val - 2.0 * e_val) / k;
131                let coil_flux = prefactor * sqrt_term * term;
132                if !coil_flux.is_finite() {
133                    return Err(FusionError::ConfigError(
134                        "vacuum field coil flux contribution became non-finite".to_string(),
135                    ));
136                }
137
138                psi_vac[[iz, ir]] += coil_flux;
139            }
140        }
141    }
142
143    if psi_vac.iter().any(|v| !v.is_finite()) {
144        return Err(FusionError::ConfigError(
145            "vacuum field result contains non-finite values".to_string(),
146        ));
147    }
148
149    Ok(psi_vac)
150}
151
152#[cfg(test)]
153mod tests {
154    use super::*;
155    use fusion_types::config::ReactorConfig;
156    use std::path::PathBuf;
157
158    fn project_root() -> PathBuf {
159        PathBuf::from(env!("CARGO_MANIFEST_DIR"))
160            .join("..")
161            .join("..")
162            .join("..")
163    }
164
165    fn config_path(relative: &str) -> String {
166        project_root().join(relative).to_string_lossy().to_string()
167    }
168
169    #[test]
170    fn test_vacuum_field_not_nan() {
171        let cfg = ReactorConfig::from_file(&config_path("iter_config.json")).unwrap();
172        let grid = cfg.create_grid();
173        let psi_vac = calculate_vacuum_field(&grid, &cfg.coils, cfg.physics.vacuum_permeability)
174            .expect("valid vacuum-field inputs");
175
176        // No NaN values
177        assert!(
178            !psi_vac.iter().any(|v| v.is_nan()),
179            "Vacuum field contains NaN"
180        );
181
182        // Should have non-zero values (coils produce flux)
183        let max_val = psi_vac.iter().cloned().fold(f64::NEG_INFINITY, f64::max);
184        assert!(max_val.abs() > 1e-10, "Vacuum field is all zeros");
185    }
186
187    #[test]
188    fn test_vacuum_field_symmetry() {
189        // Single coil at Z=0 should produce field symmetric about Z=0
190        let grid = Grid2D::new(33, 33, 1.0, 9.0, -5.0, 5.0);
191        let coil = CoilConfig {
192            name: "test".to_string(),
193            r: 5.0,
194            z: 0.0,
195            current: 1.0,
196        };
197        let psi = calculate_vacuum_field(&grid, &[coil], 1.0).expect("valid vacuum-field inputs");
198
199        // Check symmetry: psi[iz, ir] ≈ psi[nz-1-iz, ir]
200        for ir in 0..33 {
201            for iz in 0..16 {
202                let diff = (psi[[iz, ir]] - psi[[32 - iz, ir]]).abs();
203                assert!(
204                    diff < 1e-10,
205                    "Symmetry broken at iz={iz}, ir={ir}: {} vs {}",
206                    psi[[iz, ir]],
207                    psi[[32 - iz, ir]]
208                );
209            }
210        }
211    }
212
213    #[test]
214    fn test_vacuum_field_single_coil_shape() {
215        // Flux from a single coil should peak near the coil position
216        let grid = Grid2D::new(33, 33, 1.0, 9.0, -5.0, 5.0);
217        let coil = CoilConfig {
218            name: "test".to_string(),
219            r: 5.0,
220            z: 0.0,
221            current: 1.0,
222        };
223        let psi = calculate_vacuum_field(&grid, &[coil], 1.0).expect("valid vacuum-field inputs");
224
225        // Find maximum
226        let mut max_val = f64::NEG_INFINITY;
227        let mut max_iz = 0;
228        let mut max_ir = 0;
229        for iz in 0..33 {
230            for ir in 0..33 {
231                if psi[[iz, ir]] > max_val {
232                    max_val = psi[[iz, ir]];
233                    max_iz = iz;
234                    max_ir = ir;
235                }
236            }
237        }
238
239        // Maximum should be near Z=0 (middle row)
240        assert!(
241            (max_iz as f64 - 16.0).abs() < 3.0,
242            "Peak not near Z=0: iz={max_iz}"
243        );
244        // Maximum should be near R=5.0 (middle of grid)
245        assert!(
246            (max_ir as f64 - 16.0).abs() < 3.0,
247            "Peak not near R=5: ir={max_ir}"
248        );
249    }
250
251    #[test]
252    fn test_vacuum_field_rejects_invalid_runtime_inputs() {
253        let grid = Grid2D::new(17, 17, 1.0, 9.0, -5.0, 5.0);
254        let bad_coil = CoilConfig {
255            name: "bad".to_string(),
256            r: 0.0,
257            z: 0.0,
258            current: 1.0,
259        };
260        assert!(calculate_vacuum_field(&grid, &[bad_coil], 1.0).is_err());
261        assert!(calculate_vacuum_field(&grid, &[], f64::NAN).is_err());
262    }
263}
fusion_core/vacuum.rs

fusion_core/
vacuum.rs
