fusion_core/
pedestal.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 — EPED-Like Pedestal Model
7//! Simplified EPED-like pedestal model with ELM trigger logic.
8
9use fusion_types::error::{FusionError, FusionResult};
10use fusion_types::state::RadialProfiles;
11
12const MIN_T_KEV: f64 = 0.05;
13const MIN_NE: f64 = 1e-4;
14
15#[derive(Debug, Clone, Copy)]
16pub struct PedestalConfig {
17    pub beta_p_ped: f64,
18    pub rho_s: f64,
19    pub r_major: f64,
20    pub alpha_crit: f64,
21    pub tau_elm: f64,
22}
23
24impl Default for PedestalConfig {
25    fn default() -> Self {
26        Self {
27            beta_p_ped: 0.35,
28            rho_s: 2.0e-3,
29            r_major: 6.2,
30            alpha_crit: 2.5,
31            tau_elm: 1.0e-3,
32        }
33    }
34}
35
36#[derive(Debug, Clone)]
37pub struct PedestalModel {
38    pub config: PedestalConfig,
39    last_gradient: f64,
40    elm_cooldown_s: f64,
41}
42
43impl PedestalModel {
44    pub fn new(config: PedestalConfig) -> FusionResult<Self> {
45        if !config.beta_p_ped.is_finite() || config.beta_p_ped < 0.0 {
46            return Err(FusionError::ConfigError(
47                "pedestal beta_p_ped must be finite and >= 0".to_string(),
48            ));
49        }
50        if !config.rho_s.is_finite() || config.rho_s <= 0.0 {
51            return Err(FusionError::ConfigError(
52                "pedestal rho_s must be finite and > 0".to_string(),
53            ));
54        }
55        if !config.r_major.is_finite() || config.r_major <= 0.0 {
56            return Err(FusionError::ConfigError(
57                "pedestal r_major must be finite and > 0".to_string(),
58            ));
59        }
60        if !config.alpha_crit.is_finite() || config.alpha_crit <= 0.0 {
61            return Err(FusionError::ConfigError(
62                "pedestal alpha_crit must be finite and > 0".to_string(),
63            ));
64        }
65        if !config.tau_elm.is_finite() || config.tau_elm <= 0.0 {
66            return Err(FusionError::ConfigError(
67                "pedestal tau_elm must be finite and > 0".to_string(),
68            ));
69        }
70        Ok(Self {
71            config,
72            last_gradient: 0.0,
73            elm_cooldown_s: 0.0,
74        })
75    }
76
77    /// EPED-inspired width scaling: Δ_ped ~ sqrt(beta_p,ped) * (rho_s / R)
78    pub fn pedestal_width(&self) -> f64 {
79        let beta = self.config.beta_p_ped;
80        let rho_s = self.config.rho_s;
81        let r_major = self.config.r_major;
82        (beta.sqrt() * (rho_s / r_major)).clamp(0.03, 0.12)
83    }
84
85    /// Effective ELM trigger threshold: alpha_crit ~ 2.5 * s_ped.
86    fn effective_alpha_crit(&self) -> f64 {
87        let shear_proxy = 1.0 / self.pedestal_width().max(1e-4);
88        self.config.alpha_crit.max(2.5 * shear_proxy)
89    }
90
91    pub fn is_elm_triggered(&self, pressure_gradient: f64) -> bool {
92        self.elm_cooldown_s <= 0.0 && pressure_gradient.abs() > self.effective_alpha_crit()
93    }
94
95    /// Advance ELM cooldown timer.
96    pub fn advance(&mut self, dt_s: f64) {
97        self.elm_cooldown_s = (self.elm_cooldown_s - dt_s.max(0.0)).max(0.0);
98    }
99
100    pub fn last_gradient(&self) -> f64 {
101        self.last_gradient
102    }
103
104    /// Apply a fast ELM crash to pedestal region profiles.
105    pub fn apply_elm_crash(&mut self, profiles: &mut RadialProfiles) {
106        let width = self.pedestal_width();
107        let rho_start = (1.0 - width).clamp(0.75, 0.995);
108        let tau = self.config.tau_elm;
109        let burst_fraction = (1.0 - (-1.0e-3 / tau).exp()).clamp(0.05, 0.95);
110
111        for i in 0..profiles.rho.len() {
112            let rho = profiles.rho[i];
113            if rho >= rho_start {
114                let edge_w = ((rho - rho_start) / width).clamp(0.0, 1.0);
115                let temperature_drop = (0.2 + 0.6 * edge_w * burst_fraction).clamp(0.0, 0.95);
116                let density_drop = 0.5 * temperature_drop;
117
118                profiles.te[i] = (profiles.te[i] * (1.0 - temperature_drop)).max(MIN_T_KEV);
119                profiles.ti[i] = (profiles.ti[i] * (1.0 - temperature_drop)).max(MIN_T_KEV);
120                profiles.ne[i] = (profiles.ne[i] * (1.0 - density_drop)).max(MIN_NE);
121            }
122        }
123
124        self.elm_cooldown_s = 3.0 * tau;
125    }
126
127    pub fn record_gradient(&mut self, pressure_gradient: f64) {
128        self.last_gradient = pressure_gradient;
129    }
130}
131
132impl Default for PedestalModel {
133    fn default() -> Self {
134        Self::new(PedestalConfig::default()).expect("default pedestal config must be valid")
135    }
136}
137
138#[cfg(test)]
139mod tests {
140    use super::*;
141    use ndarray::Array1;
142
143    fn sample_profiles() -> RadialProfiles {
144        let n = 50;
145        let rho = Array1::linspace(0.0, 1.0, n);
146        let te = Array1::from_shape_fn(n, |i| 8.0 - 6.5 * rho[i]);
147        let ti = te.clone();
148        let ne = Array1::from_shape_fn(n, |i| 9.0 - 5.0 * rho[i]);
149        let n_impurity = Array1::zeros(n);
150        RadialProfiles {
151            rho,
152            te,
153            ti,
154            ne,
155            n_impurity,
156        }
157    }
158
159    #[test]
160    fn test_pedestal_width_eped_scaling() {
161        let model = PedestalModel::new(PedestalConfig {
162            beta_p_ped: 0.5,
163            rho_s: 3.0e-3,
164            r_major: 6.0,
165            alpha_crit: 2.5,
166            tau_elm: 1.0e-3,
167        })
168        .expect("valid pedestal config");
169        let width = model.pedestal_width();
170        assert!(width > 0.0);
171        assert!(width <= 0.12);
172    }
173
174    #[test]
175    fn test_elm_trigger_threshold() {
176        let model = PedestalModel::default();
177        assert!(!model.is_elm_triggered(1.0));
178        assert!(model.is_elm_triggered(5_000.0));
179    }
180
181    #[test]
182    fn test_apply_elm_crash_reduces_edge_profiles() {
183        let mut model = PedestalModel::default();
184        let mut profiles = sample_profiles();
185
186        let i_edge = profiles.rho.len() - 2;
187        let te_before = profiles.te[i_edge];
188        let ne_before = profiles.ne[i_edge];
189
190        model.apply_elm_crash(&mut profiles);
191
192        assert!(profiles.te[i_edge] < te_before);
193        assert!(profiles.ne[i_edge] < ne_before);
194    }
195
196    #[test]
197    fn test_elm_cooldown_blocks_retrigger_until_advanced() {
198        let mut model = PedestalModel::default();
199        let mut profiles = sample_profiles();
200
201        assert!(model.is_elm_triggered(5_000.0));
202        model.apply_elm_crash(&mut profiles);
203        assert!(
204            !model.is_elm_triggered(5_000.0),
205            "ELM cooldown should block immediate retriggering"
206        );
207
208        model.advance(10.0);
209        assert!(
210            model.is_elm_triggered(5_000.0),
211            "ELM trigger should recover after cooldown elapses"
212        );
213    }
214
215    #[test]
216    fn test_elm_crash_respects_minimum_profile_floors() {
217        let mut model = PedestalModel::default();
218        let n = 16;
219        let mut profiles = RadialProfiles {
220            rho: Array1::linspace(0.0, 1.0, n),
221            te: Array1::from_elem(n, 0.06),
222            ti: Array1::from_elem(n, 0.06),
223            ne: Array1::from_elem(n, 2.0e-4),
224            n_impurity: Array1::zeros(n),
225        };
226
227        model.apply_elm_crash(&mut profiles);
228        assert!(profiles.te.iter().all(|v| *v >= MIN_T_KEV));
229        assert!(profiles.ti.iter().all(|v| *v >= MIN_T_KEV));
230        assert!(profiles.ne.iter().all(|v| *v >= MIN_NE));
231    }
232
233    #[test]
234    fn test_pedestal_constructor_rejects_invalid_config() {
235        let bad = [
236            PedestalConfig {
237                beta_p_ped: -1.0,
238                ..Default::default()
239            },
240            PedestalConfig {
241                rho_s: 0.0,
242                ..Default::default()
243            },
244            PedestalConfig {
245                r_major: 0.0,
246                ..Default::default()
247            },
248            PedestalConfig {
249                alpha_crit: 0.0,
250                ..Default::default()
251            },
252            PedestalConfig {
253                tau_elm: 0.0,
254                ..Default::default()
255            },
256        ];
257        for cfg in bad {
258            let err = PedestalModel::new(cfg).expect_err("invalid config must error");
259            match err {
260                FusionError::ConfigError(msg) => {
261                    assert!(msg.contains("pedestal"));
262                }
263                other => panic!("Unexpected error: {other:?}"),
264            }
265        }
266    }
267}
fusion_core/pedestal.rs

fusion_core/
pedestal.rs
