#![warn(missing_docs)]
pub mod alara;
pub mod armi;
pub mod fluka;
pub mod mcnp;
pub mod partisn;
pub mod serpent;
use nucleide_material::Material;
use nucleide_nuclei::NuclideId;
#[derive(Debug, Clone, Copy, PartialEq, Eq, Hash)]
pub enum Code {
Mcnp,
Serpent,
Fluka,
Alara,
Partisn,
}
impl Code {
pub fn all() -> [Code; 5] {
[
Code::Mcnp,
Code::Serpent,
Code::Fluka,
Code::Alara,
Code::Partisn,
]
}
}
impl std::fmt::Display for Code {
fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result {
match self {
Code::Mcnp => write!(f, "MCNP"),
Code::Serpent => write!(f, "Serpent"),
Code::Fluka => write!(f, "FLUKA"),
Code::Alara => write!(f, "ALARA"),
Code::Partisn => write!(f, "PARTISN"),
}
}
}
#[derive(Debug, Clone)]
pub struct EmitOptions {
pub name: String,
pub mcnp_number: u32,
pub xs_suffix: String,
pub density: Option<f64>,
pub serpent_lib: String,
pub fluka_fid: u32,
pub partisn_zone: u32,
}
impl EmitOptions {
pub fn new(name: impl Into<String>) -> Self {
Self {
name: name.into(),
mcnp_number: 1,
xs_suffix: "80c".to_string(),
density: None,
serpent_lib: "03c".to_string(),
fluka_fid: 1,
partisn_zone: 1,
}
}
pub fn with_density(mut self, density: f64) -> Self {
self.density = Some(density);
self
}
fn density_for(&self, mat: &Material) -> Result<f64> {
self.density
.or_else(|| mat.density())
.ok_or(Error::MissingDensity)
}
}
#[derive(Debug, Clone, PartialEq)]
pub struct Dropped {
pub id: NuclideId,
pub mass: f64,
pub reason: String,
}
#[derive(Debug, Clone)]
pub struct Emitted {
pub code: Code,
pub text: String,
pub accounted: Vec<(NuclideId, f64)>,
pub dropped: Vec<Dropped>,
pub reparsed: bool,
}
impl Emitted {
pub fn mass_out(&self) -> f64 {
self.accounted.iter().map(|(_, m)| m).sum()
}
}
#[derive(Debug, Clone)]
pub struct DriftRow {
pub code: Code,
pub mass_in: f64,
pub mass_out: f64,
pub rel_drift: f64,
pub dropped: Vec<Dropped>,
pub reparsed: bool,
}
#[derive(Debug, Clone)]
pub struct DriftTable {
pub name: String,
pub rows: Vec<DriftRow>,
}
impl DriftTable {
pub fn worst_rel_drift(&self) -> f64 {
self.rows
.iter()
.map(|r| r.rel_drift.abs())
.fold(0.0, f64::max)
}
}
#[derive(Debug, thiserror::Error)]
#[non_exhaustive]
pub enum Error {
#[error("material is empty or its masses sum to a non-positive value")]
Degenerate,
#[error("emission requires a mass density but none was set")]
MissingDensity,
#[error("emitted {code} text failed to re-parse: {detail}")]
Reparse {
code: Code,
detail: String,
},
#[error(transparent)]
Material(#[from] nucleide_material::Error),
#[error(transparent)]
Nuclei(#[from] nucleide_nuclei::Error),
#[error(transparent)]
Mcnp(#[from] nucleide_mcnp_io::inp::Error),
#[error(transparent)]
Alara(#[from] nucleide_alara_io::Error),
#[error(transparent)]
Fluka(#[from] nucleide_fluka_io::material::Error),
#[error("invalid ARMI key `{key}`: {reason}")]
ArmiKey {
key: String,
reason: String,
},
}
pub type Result<T> = std::result::Result<T, Error>;
pub fn emit_all(mat: &Material, opts: &EmitOptions) -> Result<Vec<Emitted>> {
if mat.mass() <= 0.0 || !mat.mass().is_finite() {
return Err(Error::Degenerate);
}
Ok(vec![
mcnp::emit_mcnp(mat, opts)?,
serpent::emit_serpent(mat, opts)?,
fluka::emit_fluka(mat, opts)?,
alara::emit_alara(mat, opts)?,
partisn::emit_partisn(mat, opts)?,
])
}
pub fn drift_table(mat: &Material, emitted: &[Emitted]) -> DriftTable {
let mass_in = mat.mass();
let rows = emitted
.iter()
.map(|e| {
let mass_out = e.mass_out();
DriftRow {
code: e.code,
mass_in,
mass_out,
rel_drift: if mass_in == 0.0 {
0.0
} else {
(mass_in - mass_out) / mass_in
},
dropped: e.dropped.clone(),
reparsed: e.reparsed,
}
})
.collect();
DriftTable {
name: String::new(),
rows,
}
}
pub fn emit_drift(mat: &Material, opts: &EmitOptions) -> Result<(Vec<Emitted>, DriftTable)> {
let emitted = emit_all(mat, opts)?;
let mut table = drift_table(mat, &emitted);
table.name = opts.name.clone();
Ok((emitted, table))
}
#[cfg(test)]
mod tests {
use super::*;
fn metal() -> Material {
let mut mat = Material::new();
mat.add_nuclide(NuclideId::from_name("U235").unwrap(), 5.0);
mat.add_nuclide(NuclideId::from_name("U238").unwrap(), 95.0);
mat
}
#[test]
fn emit_all_covers_five_codes_lossless() {
let opts = EmitOptions::new("umetal").with_density(19.1);
let (emitted, table) = emit_drift(&metal(), &opts).unwrap();
assert_eq!(emitted.len(), 5);
assert_eq!(
emitted.iter().map(|e| e.code).collect::<Vec<_>>(),
Code::all()
);
assert_eq!(table.name, "umetal");
assert_eq!(table.rows.len(), 5);
for row in &table.rows {
assert!((row.mass_in - 100.0).abs() < 1e-12, "{:?}", row.code);
assert!((row.mass_out - 100.0).abs() < 1e-9, "{:?}", row.code);
assert!(row.rel_drift.abs() < 1e-9, "{:?}", row.code);
assert!(row.dropped.is_empty(), "{:?}", row.code);
}
assert!(table.worst_rel_drift() < 1e-9);
let reparsed: Vec<Code> = emitted
.iter()
.filter(|e| e.reparsed)
.map(|e| e.code)
.collect();
assert_eq!(reparsed, vec![Code::Mcnp, Code::Alara]);
}
#[test]
fn drift_reports_fluka_loss() {
let mut mat = Material::new();
mat.add_nuclide(NuclideId::from_name("H1").unwrap(), 20.0);
mat.add_nuclide(NuclideId::from_name("O16").unwrap(), 80.0);
let opts = EmitOptions::new("water").with_density(1.0);
let (_, table) = emit_drift(&mat, &opts).unwrap();
let fluka = table.rows.iter().find(|r| r.code == Code::Fluka).unwrap();
assert!((fluka.rel_drift - 0.8).abs() < 1e-12);
assert_eq!(fluka.dropped.len(), 1);
assert_eq!(fluka.dropped[0].id, NuclideId::from_name("O16").unwrap());
assert!((table.worst_rel_drift() - 0.8).abs() < 1e-12);
for row in table.rows.iter().filter(|r| r.code != Code::Fluka) {
assert!(row.rel_drift.abs() < 1e-9, "{:?}", row.code);
}
}
#[test]
fn code_display_names() {
assert_eq!(
Code::all()
.iter()
.map(|c| c.to_string())
.collect::<Vec<_>>(),
vec!["MCNP", "Serpent", "FLUKA", "ALARA", "PARTISN"]
);
}
#[test]
fn empty_material_is_degenerate() {
let opts = EmitOptions::new("void").with_density(1.0);
assert!(matches!(
emit_all(&Material::new(), &opts),
Err(Error::Degenerate)
));
}
}