Skip to main content

nucleide_emit/
lib.rs

1#![warn(missing_docs)]
2//! Single-material emission to legacy transport-code cards.
3//!
4//! [`emit_all`] renders one [`Material`] through five code dialects — MCNP,
5//! Serpent, FLUKA, ALARA, PARTISN — and [`drift_table`] reports how much mass
6//! survives each translation. Every emitter is pure glue over the workspace
7//! `*-io` crates plus [`nucleide_nuclei::dialects`]; no new physics lives here.
8//!
9//! # Drift semantics
10//!
11//! Each emitter returns the mass it could represent per nuclide
12//! ([`Emitted::accounted`]) plus what it had to drop ([`Emitted::dropped`],
13//! e.g. a nuclide with no FLUKA name). [`DriftRow::rel_drift`] is
14//! `(mass_in - mass_out) / mass_in`. [`Emitted::reparsed`] tells whether the
15//! emitted text was machine-verified by feeding it back through that code's
16//! reader (MCNP and ALARA only — Serpent, FLUKA, and PARTISN have no material
17//! readers in this workspace, so their drift is analytic).
18//!
19//! # Example
20//!
21//! ```rust
22//! use nucleide_emit::{EmitOptions, emit_drift};
23//! use nucleide_material::Material;
24//! use nucleide_nuclei::NuclideId;
25//!
26//! let mut mat = Material::new();
27//! mat.add_nuclide(NuclideId::from_name("U235").unwrap(), 5.0);
28//! mat.add_nuclide(NuclideId::from_name("U238").unwrap(), 95.0);
29//! let opts = EmitOptions::new("leu").with_density(10.0);
30//! let (emitted, table) = emit_drift(&mat, &opts).unwrap();
31//! assert_eq!(emitted.len(), 5);
32//! assert!(table.worst_rel_drift() < 1e-9);
33//! ```
34
35pub mod alara;
36pub mod armi;
37pub mod fluka;
38pub mod mcnp;
39pub mod partisn;
40pub mod serpent;
41
42use nucleide_material::Material;
43use nucleide_nuclei::NuclideId;
44
45/// One of the five supported transport-code dialects.
46#[derive(Debug, Clone, Copy, PartialEq, Eq, Hash)]
47pub enum Code {
48    /// MCNP `m<number>` material card.
49    Mcnp,
50    /// Serpent `mat <name> <-density>` material card.
51    Serpent,
52    /// FLUKA `COMPOUND` card.
53    Fluka,
54    /// ALARA `mixture` block.
55    Alara,
56    /// PARTISN single-zone deck.
57    Partisn,
58}
59
60impl Code {
61    /// All codes in emission order.
62    pub fn all() -> [Code; 5] {
63        [
64            Code::Mcnp,
65            Code::Serpent,
66            Code::Fluka,
67            Code::Alara,
68            Code::Partisn,
69        ]
70    }
71}
72
73impl std::fmt::Display for Code {
74    fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result {
75        match self {
76            Code::Mcnp => write!(f, "MCNP"),
77            Code::Serpent => write!(f, "Serpent"),
78            Code::Fluka => write!(f, "FLUKA"),
79            Code::Alara => write!(f, "ALARA"),
80            Code::Partisn => write!(f, "PARTISN"),
81        }
82    }
83}
84
85/// Shared emission settings.
86#[derive(Debug, Clone)]
87pub struct EmitOptions {
88    /// Material/card name used by every dialect.
89    pub name: String,
90    /// MCNP `m<number>` card number.
91    pub mcnp_number: u32,
92    /// MCNP cross-section library suffix (`80c` renders `92235.80c`).
93    pub xs_suffix: String,
94    /// Mass density [g/cm³] for dialects that need one (Serpent, FLUKA,
95    /// PARTISN). Falls back to [`Material::density`]; errors when absent.
96    pub density: Option<f64>,
97    /// Serpent cross-section library suffix (`03c` renders `92235.03c`).
98    pub serpent_lib: String,
99    /// FLUKA material index number.
100    pub fluka_fid: u32,
101    /// PARTISN zone id for the single emitted zone.
102    pub partisn_zone: u32,
103}
104
105impl EmitOptions {
106    /// Options with conventional defaults (`m1`, `80c`, FLUKA fid 1, zone 1).
107    pub fn new(name: impl Into<String>) -> Self {
108        Self {
109            name: name.into(),
110            mcnp_number: 1,
111            xs_suffix: "80c".to_string(),
112            density: None,
113            serpent_lib: "03c".to_string(),
114            fluka_fid: 1,
115            partisn_zone: 1,
116        }
117    }
118
119    /// Override the mass density [g/cm³].
120    pub fn with_density(mut self, density: f64) -> Self {
121        self.density = Some(density);
122        self
123    }
124
125    fn density_for(&self, mat: &Material) -> Result<f64> {
126        self.density
127            .or_else(|| mat.density())
128            .ok_or(Error::MissingDensity)
129    }
130}
131
132/// One nuclide the emitter could not represent.
133#[derive(Debug, Clone, PartialEq)]
134pub struct Dropped {
135    /// Nuclide that was skipped.
136    pub id: NuclideId,
137    /// Mass \[g\] left out of the emitted cards.
138    pub mass: f64,
139    /// Machine-readable reason (e.g. `"no-fluka-name"`, `"no-atomic-mass"`).
140    pub reason: String,
141}
142
143/// Cards emitted for one dialect.
144#[derive(Debug, Clone)]
145pub struct Emitted {
146    /// Dialect these cards belong to.
147    pub code: Code,
148    /// Card text, ready to paste into a deck.
149    pub text: String,
150    /// Mass \[g\] accounted per emitted nuclide, in [`Material`] order.
151    pub accounted: Vec<(NuclideId, f64)>,
152    /// Nuclides skipped with reasons.
153    pub dropped: Vec<Dropped>,
154    /// Whether `text` was verified by re-parsing with the code's own reader.
155    /// Only MCNP and ALARA have material readers in this workspace.
156    pub reparsed: bool,
157}
158
159impl Emitted {
160    /// Mass \[g\] represented by these cards.
161    pub fn mass_out(&self) -> f64 {
162        self.accounted.iter().map(|(_, m)| m).sum()
163    }
164}
165
166/// One row of the mass-drift report.
167#[derive(Debug, Clone)]
168pub struct DriftRow {
169    /// Dialect this row covers.
170    pub code: Code,
171    /// Input mass \[g\].
172    pub mass_in: f64,
173    /// Mass represented in the emitted cards \[g\].
174    pub mass_out: f64,
175    /// `(mass_in - mass_out) / mass_in`; zero for lossless emission.
176    pub rel_drift: f64,
177    /// Nuclides skipped with reasons.
178    pub dropped: Vec<Dropped>,
179    /// Whether the emitted text was re-parse verified.
180    pub reparsed: bool,
181}
182
183/// Mass conservation across all five dialects.
184#[derive(Debug, Clone)]
185pub struct DriftTable {
186    /// Material name from [`EmitOptions::name`].
187    pub name: String,
188    /// One row per dialect, in [`Code::all`] order.
189    pub rows: Vec<DriftRow>,
190}
191
192impl DriftTable {
193    /// Largest relative drift over all dialects.
194    pub fn worst_rel_drift(&self) -> f64 {
195        self.rows
196            .iter()
197            .map(|r| r.rel_drift.abs())
198            .fold(0.0, f64::max)
199    }
200}
201
202/// Errors raised while emitting cards.
203#[derive(Debug, thiserror::Error)]
204#[non_exhaustive]
205pub enum Error {
206    /// Material is empty or its masses sum to a non-positive value.
207    #[error("material is empty or its masses sum to a non-positive value")]
208    Degenerate,
209    /// A dialect needs a mass density but neither options nor material has one.
210    #[error("emission requires a mass density but none was set")]
211    MissingDensity,
212    /// The emitted text failed to re-parse with the code's own reader.
213    #[error("emitted {code} text failed to re-parse: {detail}")]
214    Reparse {
215        /// Dialect whose text did not survive its reader.
216        code: Code,
217        /// Reader error.
218        detail: String,
219    },
220    /// Material-layer failure (mass tables, fractions).
221    #[error(transparent)]
222    Material(#[from] nucleide_material::Error),
223    /// Nuclide-dialect failure.
224    #[error(transparent)]
225    Nuclei(#[from] nucleide_nuclei::Error),
226    /// MCNP reader failure during re-parse verification.
227    #[error(transparent)]
228    Mcnp(#[from] nucleide_mcnp_io::inp::Error),
229    /// ALARA reader failure during re-parse verification.
230    #[error(transparent)]
231    Alara(#[from] nucleide_alara_io::Error),
232    /// FLUKA naming failure.
233    #[error(transparent)]
234    Fluka(#[from] nucleide_fluka_io::material::Error),
235    /// An ARMI mass-fraction key (or its value) was rejected by the v1
236    /// caller-side rules in [`armi`]: elemental keys, bare `AM242`, and
237    /// negative/non-finite masses are caller errors.
238    #[error("invalid ARMI key `{key}`: {reason}")]
239    ArmiKey {
240        /// The offending ARMI-side key.
241        key: String,
242        /// Why it was rejected (expand-first, disambiguation, value range).
243        reason: String,
244    },
245}
246
247/// Crate-local result alias.
248pub type Result<T> = std::result::Result<T, Error>;
249
250/// Emit one [`Material`] through all five dialects, in [`Code::all`] order.
251pub fn emit_all(mat: &Material, opts: &EmitOptions) -> Result<Vec<Emitted>> {
252    if mat.mass() <= 0.0 || !mat.mass().is_finite() {
253        return Err(Error::Degenerate);
254    }
255    Ok(vec![
256        mcnp::emit_mcnp(mat, opts)?,
257        serpent::emit_serpent(mat, opts)?,
258        fluka::emit_fluka(mat, opts)?,
259        alara::emit_alara(mat, opts)?,
260        partisn::emit_partisn(mat, opts)?,
261    ])
262}
263
264/// Build the mass-drift report for already-emitted cards.
265pub fn drift_table(mat: &Material, emitted: &[Emitted]) -> DriftTable {
266    let mass_in = mat.mass();
267    let rows = emitted
268        .iter()
269        .map(|e| {
270            let mass_out = e.mass_out();
271            DriftRow {
272                code: e.code,
273                mass_in,
274                mass_out,
275                rel_drift: if mass_in == 0.0 {
276                    0.0
277                } else {
278                    (mass_in - mass_out) / mass_in
279                },
280                dropped: e.dropped.clone(),
281                reparsed: e.reparsed,
282            }
283        })
284        .collect();
285    DriftTable {
286        name: String::new(),
287        rows,
288    }
289}
290
291/// Emit through all dialects and report mass drift in one call.
292pub fn emit_drift(mat: &Material, opts: &EmitOptions) -> Result<(Vec<Emitted>, DriftTable)> {
293    let emitted = emit_all(mat, opts)?;
294    let mut table = drift_table(mat, &emitted);
295    table.name = opts.name.clone();
296    Ok((emitted, table))
297}
298
299#[cfg(test)]
300mod tests {
301    use super::*;
302
303    fn metal() -> Material {
304        // Uranium metal: every dialect represents both nuclides, so the
305        // full table is lossless. (Light isotopes such as H1/O16 have no
306        // FLUKA isotope-table entry and exercise the dropped path instead.)
307        let mut mat = Material::new();
308        mat.add_nuclide(NuclideId::from_name("U235").unwrap(), 5.0);
309        mat.add_nuclide(NuclideId::from_name("U238").unwrap(), 95.0);
310        mat
311    }
312
313    #[test]
314    fn emit_all_covers_five_codes_lossless() {
315        let opts = EmitOptions::new("umetal").with_density(19.1);
316        let (emitted, table) = emit_drift(&metal(), &opts).unwrap();
317        assert_eq!(emitted.len(), 5);
318        assert_eq!(
319            emitted.iter().map(|e| e.code).collect::<Vec<_>>(),
320            Code::all()
321        );
322        assert_eq!(table.name, "umetal");
323        assert_eq!(table.rows.len(), 5);
324        for row in &table.rows {
325            assert!((row.mass_in - 100.0).abs() < 1e-12, "{:?}", row.code);
326            assert!((row.mass_out - 100.0).abs() < 1e-9, "{:?}", row.code);
327            assert!(row.rel_drift.abs() < 1e-9, "{:?}", row.code);
328            assert!(row.dropped.is_empty(), "{:?}", row.code);
329        }
330        assert!(table.worst_rel_drift() < 1e-9);
331        let reparsed: Vec<Code> = emitted
332            .iter()
333            .filter(|e| e.reparsed)
334            .map(|e| e.code)
335            .collect();
336        assert_eq!(reparsed, vec![Code::Mcnp, Code::Alara]);
337    }
338
339    #[test]
340    fn drift_reports_fluka_loss() {
341        // O16 has no isotope entry in the vendored FLUKA table (H1 maps to
342        // HYDROG-1), so 80 of 100 g drift away on that row while the other
343        // four stay lossless.
344        let mut mat = Material::new();
345        mat.add_nuclide(NuclideId::from_name("H1").unwrap(), 20.0);
346        mat.add_nuclide(NuclideId::from_name("O16").unwrap(), 80.0);
347        let opts = EmitOptions::new("water").with_density(1.0);
348        let (_, table) = emit_drift(&mat, &opts).unwrap();
349        let fluka = table.rows.iter().find(|r| r.code == Code::Fluka).unwrap();
350        assert!((fluka.rel_drift - 0.8).abs() < 1e-12);
351        assert_eq!(fluka.dropped.len(), 1);
352        assert_eq!(fluka.dropped[0].id, NuclideId::from_name("O16").unwrap());
353        assert!((table.worst_rel_drift() - 0.8).abs() < 1e-12);
354        for row in table.rows.iter().filter(|r| r.code != Code::Fluka) {
355            assert!(row.rel_drift.abs() < 1e-9, "{:?}", row.code);
356        }
357    }
358
359    #[test]
360    fn code_display_names() {
361        assert_eq!(
362            Code::all()
363                .iter()
364                .map(|c| c.to_string())
365                .collect::<Vec<_>>(),
366            vec!["MCNP", "Serpent", "FLUKA", "ALARA", "PARTISN"]
367        );
368    }
369
370    #[test]
371    fn empty_material_is_degenerate() {
372        let opts = EmitOptions::new("void").with_density(1.0);
373        assert!(matches!(
374            emit_all(&Material::new(), &opts),
375            Err(Error::Degenerate)
376        ));
377    }
378}