fusion_core/
ignition.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 — Ignition
7//! Ignition physics: Bosch-Hale D-T reaction rate and thermodynamics.
8//!
9//! Port of `fusion_ignition_sim.py` lines 14-107.
10//! Calculates fusion power, alpha heating, and Q-factor.
11
12use crate::kernel::FusionKernel;
13use fusion_types::constants::E_FUSION_DT;
14use fusion_types::error::{FusionError, FusionResult};
15use fusion_types::state::ThermodynamicsResult;
16
17/// Peak density [m⁻³]. Python line 56.
18const N_PEAK: f64 = 1.0e20;
19
20/// Peak temperature [keV]. Python line 57.
21const T_PEAK_KEV: f64 = 20.0;
22
23/// Density profile exponent. Python line 63.
24const DENSITY_EXPONENT: f64 = 0.5;
25
26/// Temperature profile exponent. Python line 64.
27const TEMPERATURE_EXPONENT: f64 = 1.0;
28
29/// Alpha particle energy fraction (3.5/17.6 MeV ≈ 0.199).
30/// Python line 82: `P_alpha = P_fusion * 0.2`
31const ALPHA_FRACTION: f64 = 0.2;
32
33/// Energy confinement time [s]. Python line 88.
34const TAU_E: f64 = 3.0;
35
36/// keV to Joules conversion: 1 keV = 1000 eV × 1.602e-19 J/eV
37const KEV_TO_JOULES: f64 = 1.602e-16;
38
39/// Bosch-Hale D-T fusion reaction rate ⟨σv⟩ in m³/s.
40///
41/// Uses NRL Plasma Formulary approximation:
42///   σv = 3.68e-18 / T^(2/3) × exp(-19.94 / T^(1/3))
43///
44/// Valid for T < 100 keV. Caller must provide finite T > 0.
45pub fn bosch_hale_dt(t_kev: f64) -> FusionResult<f64> {
46    if !t_kev.is_finite() || t_kev <= 0.0 {
47        return Err(FusionError::ConfigError(
48            "ignition temperature must be finite and > 0 keV".to_string(),
49        ));
50    }
51    let rate = 3.68e-18 / t_kev.powf(2.0 / 3.0) * (-19.94 / t_kev.powf(1.0 / 3.0)).exp();
52    if !rate.is_finite() {
53        return Err(FusionError::ConfigError(
54            "ignition Bosch-Hale rate became non-finite".to_string(),
55        ));
56    }
57    Ok(rate)
58}
59
60/// Calculate full thermodynamics from equilibrium state.
61///
62/// Algorithm:
63/// 1. Normalize flux: ψ_norm = (Ψ - Ψ_axis) / (Ψ_boundary - Ψ_axis)
64/// 2. Profiles: n(ψ) = n_peak·(1 - ψ²)^0.5, T(ψ) = T_peak·(1 - ψ²)^1.0
65/// 3. Fusion power: P_fus = ∫ nD·nT·⟨σv⟩·E_fus dV, dV = dR·dZ·2πR
66/// 4. Alpha heating: P_alpha = 0.2·P_fusion
67/// 5. Thermal energy: W_th = ∫ 3·n·T dV
68/// 6. Losses: P_loss = W_th / τ_E
69/// 7. Q = P_fusion / P_aux
70pub fn calculate_thermodynamics(
71    kernel: &FusionKernel,
72    p_aux_mw: f64,
73) -> FusionResult<ThermodynamicsResult> {
74    if !p_aux_mw.is_finite() || p_aux_mw < 0.0 {
75        return Err(FusionError::ConfigError(
76            "ignition p_aux_mw must be finite and >= 0".to_string(),
77        ));
78    }
79
80    let grid = kernel.grid();
81    let psi = kernel.psi();
82    let state = kernel.state();
83    let nz = grid.nz;
84    let nr = grid.nr;
85    let dr = grid.dr;
86    let dz = grid.dz;
87
88    let psi_axis = state.psi_axis;
89    let psi_boundary = state.psi_boundary;
90    let denom = psi_boundary - psi_axis;
91    if !denom.is_finite() || denom.abs() < 1e-9 {
92        return Err(FusionError::ConfigError(
93            "ignition flux normalization denominator must be finite and non-zero".to_string(),
94        ));
95    }
96
97    let mut p_fusion_w = 0.0_f64;
98    let mut w_thermal_j = 0.0_f64;
99    let mut t_peak_actual = 0.0_f64;
100
101    for iz in 0..nz {
102        for ir in 0..nr {
103            let psi_norm = (psi[[iz, ir]] - psi_axis) / denom;
104
105            if (0.0..1.0).contains(&psi_norm) {
106                let psi_norm_sq = psi_norm * psi_norm;
107                let one_minus = (1.0 - psi_norm_sq).max(0.0);
108
109                // Profiles
110                let n_e = N_PEAK * one_minus.powf(DENSITY_EXPONENT); // m⁻³
111                let t_kev = T_PEAK_KEV * one_minus.powf(TEMPERATURE_EXPONENT); // keV
112
113                if t_kev > t_peak_actual {
114                    t_peak_actual = t_kev;
115                }
116
117                // D-T: assume 50/50 mix → nD = nT = n_e / 2
118                let n_d = n_e / 2.0;
119                let n_t = n_e / 2.0;
120
121                let sigma_v = bosch_hale_dt(t_kev)?;
122
123                let r = grid.rr[[iz, ir]];
124                let dv = dr * dz * 2.0 * std::f64::consts::PI * r; // toroidal volume element
125
126                // Fusion power density × volume
127                p_fusion_w += n_d * n_t * sigma_v * E_FUSION_DT * dv;
128
129                // Thermal energy: W = 3 n T (using keV→J conversion)
130                w_thermal_j += 3.0 * n_e * t_kev * KEV_TO_JOULES * dv;
131            }
132        }
133    }
134
135    let p_fusion_mw = p_fusion_w / 1e6;
136    let p_alpha_mw = p_fusion_mw * ALPHA_FRACTION;
137    let w_thermal_mj = w_thermal_j / 1e6;
138    let p_loss_mw = w_thermal_mj / TAU_E; // MW = MJ/s
139
140    let q_factor = if p_aux_mw > 1e-9 {
141        p_fusion_mw / p_aux_mw
142    } else {
143        0.0
144    };
145
146    let net_mw = p_alpha_mw + p_aux_mw - p_loss_mw;
147
148    Ok(ThermodynamicsResult {
149        p_fusion_mw,
150        p_alpha_mw,
151        p_loss_mw,
152        p_aux_mw,
153        net_mw,
154        q_factor,
155        t_peak_kev: t_peak_actual,
156        w_thermal_mj,
157    })
158}
159
160#[cfg(test)]
161mod tests {
162    use super::*;
163
164    #[test]
165    fn test_bosch_hale_zero_temp() {
166        // At T=0.1 keV, rate should be extremely small
167        assert!(
168            bosch_hale_dt(0.1).expect("valid keV input") < 1e-30,
169            "Rate at T=0.1 keV should be tiny"
170        );
171    }
172
173    #[test]
174    fn test_bosch_hale_20kev() {
175        let rate = bosch_hale_dt(20.0).expect("valid keV input");
176        assert!(
177            rate > 1e-22 && rate < 1e-21,
178            "Rate at 20keV: {rate}, expected O(1e-22)"
179        );
180    }
181
182    #[test]
183    fn test_bosch_hale_peak_higher() {
184        let rate_60 = bosch_hale_dt(60.0).expect("valid keV input");
185        let rate_1 = bosch_hale_dt(1.0).expect("valid keV input");
186        assert!(
187            rate_60 > rate_1 * 1000.0,
188            "Rate at 60keV ({rate_60}) should be >> rate at 1keV ({rate_1})"
189        );
190    }
191
192    #[test]
193    fn test_bosch_hale_monotonic_rise() {
194        // Rate should increase monotonically from 1 to ~60 keV
195        let mut prev = bosch_hale_dt(1.0).expect("valid keV input");
196        for t in [5.0, 10.0, 20.0, 40.0, 60.0] {
197            let rate = bosch_hale_dt(t).expect("valid keV input");
198            assert!(
199                rate > prev,
200                "Rate should increase: T={t}, rate={rate}, prev={prev}"
201            );
202            prev = rate;
203        }
204    }
205
206    #[test]
207    fn test_thermodynamics_after_solve() {
208        use std::path::PathBuf;
209
210        fn project_root() -> PathBuf {
211            PathBuf::from(env!("CARGO_MANIFEST_DIR"))
212                .join("..")
213                .join("..")
214                .join("..")
215        }
216        fn config_path(relative: &str) -> String {
217            project_root().join(relative).to_string_lossy().to_string()
218        }
219
220        // Use smaller grid for faster test
221        let mut kernel =
222            FusionKernel::from_file(&config_path("validation/iter_validated_config.json")).unwrap();
223        kernel.solve_equilibrium().unwrap();
224
225        let result = calculate_thermodynamics(&kernel, 50.0).expect("valid thermodynamics input");
226
227        // Fusion power should be positive and finite
228        assert!(
229            result.p_fusion_mw > 0.0,
230            "P_fusion = {}",
231            result.p_fusion_mw
232        );
233        assert!(result.p_fusion_mw.is_finite(), "P_fusion must be finite");
234
235        // Q-factor should be positive
236        assert!(result.q_factor > 0.0, "Q = {}", result.q_factor);
237
238        // Alpha power should be 20% of fusion
239        let alpha_ratio = result.p_alpha_mw / result.p_fusion_mw;
240        assert!(
241            (alpha_ratio - 0.2).abs() < 1e-10,
242            "Alpha ratio = {alpha_ratio}, expected 0.2"
243        );
244    }
245
246    #[test]
247    fn test_ignition_rejects_invalid_inputs() {
248        assert!(bosch_hale_dt(0.0).is_err());
249        assert!(bosch_hale_dt(f64::NAN).is_err());
250
251        use std::path::PathBuf;
252        fn project_root() -> PathBuf {
253            PathBuf::from(env!("CARGO_MANIFEST_DIR"))
254                .join("..")
255                .join("..")
256                .join("..")
257        }
258        fn config_path(relative: &str) -> String {
259            project_root().join(relative).to_string_lossy().to_string()
260        }
261
262        let kernel =
263            FusionKernel::from_file(&config_path("validation/iter_validated_config.json")).unwrap();
264        assert!(calculate_thermodynamics(&kernel, f64::NAN).is_err());
265    }
266}
fusion_core/ignition.rs

fusion_core/
ignition.rs
