fusion_core/
pedestal.rs

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

fusion_core/
pedestal.rs
