1#![warn(missing_docs)]
2pub 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#[derive(Debug, Clone, Copy, PartialEq, Eq, Hash)]
47pub enum Code {
48 Mcnp,
50 Serpent,
52 Fluka,
54 Alara,
56 Partisn,
58}
59
60impl Code {
61 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#[derive(Debug, Clone)]
87pub struct EmitOptions {
88 pub name: String,
90 pub mcnp_number: u32,
92 pub xs_suffix: String,
94 pub density: Option<f64>,
97 pub serpent_lib: String,
99 pub fluka_fid: u32,
101 pub partisn_zone: u32,
103}
104
105impl EmitOptions {
106 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 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#[derive(Debug, Clone, PartialEq)]
134pub struct Dropped {
135 pub id: NuclideId,
137 pub mass: f64,
139 pub reason: String,
141}
142
143#[derive(Debug, Clone)]
145pub struct Emitted {
146 pub code: Code,
148 pub text: String,
150 pub accounted: Vec<(NuclideId, f64)>,
152 pub dropped: Vec<Dropped>,
154 pub reparsed: bool,
157}
158
159impl Emitted {
160 pub fn mass_out(&self) -> f64 {
162 self.accounted.iter().map(|(_, m)| m).sum()
163 }
164}
165
166#[derive(Debug, Clone)]
168pub struct DriftRow {
169 pub code: Code,
171 pub mass_in: f64,
173 pub mass_out: f64,
175 pub rel_drift: f64,
177 pub dropped: Vec<Dropped>,
179 pub reparsed: bool,
181}
182
183#[derive(Debug, Clone)]
185pub struct DriftTable {
186 pub name: String,
188 pub rows: Vec<DriftRow>,
190}
191
192impl DriftTable {
193 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#[derive(Debug, thiserror::Error)]
204#[non_exhaustive]
205pub enum Error {
206 #[error("material is empty or its masses sum to a non-positive value")]
208 Degenerate,
209 #[error("emission requires a mass density but none was set")]
211 MissingDensity,
212 #[error("emitted {code} text failed to re-parse: {detail}")]
214 Reparse {
215 code: Code,
217 detail: String,
219 },
220 #[error(transparent)]
222 Material(#[from] nucleide_material::Error),
223 #[error(transparent)]
225 Nuclei(#[from] nucleide_nuclei::Error),
226 #[error(transparent)]
228 Mcnp(#[from] nucleide_mcnp_io::inp::Error),
229 #[error(transparent)]
231 Alara(#[from] nucleide_alara_io::Error),
232 #[error(transparent)]
234 Fluka(#[from] nucleide_fluka_io::material::Error),
235 #[error("invalid ARMI key `{key}`: {reason}")]
239 ArmiKey {
240 key: String,
242 reason: String,
244 },
245}
246
247pub type Result<T> = std::result::Result<T, Error>;
249
250pub 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
264pub 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
291pub 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 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 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}