Survey ways to define multi-element materials

import openmc
import matplotlib.pyplot as plt

from openmc_workbench_config import DATA_ROOT

xs_xml = DATA_ROOT / 'endf-b8.1' / 'cross_sections.xml'
openmc.config['cross_sections'] = xs_xml


def plotxs(dic):
    fig, axes = plt.subplots(1, 3, figsize=(15, 5))

    for ax, (name, mat) in zip(axes, dic.items()):
        openmc.plot_xs(
            reactions={mat: ['total']},
            axis=ax,
        )
        ax.set_title(name)
        ax.get_legend().remove()
        ax.lines[-1].set(color='g', marker='.', ls='None', ms=1)

    fig.tight_layout()

steel = openmc.Material(name='Steel, Stainless 316L')
steel.add_element('C', 0.001384)
steel.add_element('Mn', 0.020165)
steel.add_element('P', 0.000805)
steel.add_element('S', 0.000518)
steel.add_element('Si', 0.019722)
steel.add_element('Cr', 0.181098)
steel.add_element('Ni', 0.113247)
steel.add_element('Mo', 0.014432)
steel.add_element('Fe', 0.648629)
steel.set_density('g/cm3', 8.0)

concrete = openmc.Material(name='Concrete, Ordinary (NIST)')
concrete.add_element('H', 0.305322)
concrete.add_element('C', 0.002880)
concrete.add_element('O', 0.500411)
concrete.add_element('Na', 0.009212)
concrete.add_element('Mg', 0.000725)
concrete.add_element('Al', 0.010298)
concrete.add_element('Si', 0.151046)
concrete.add_element('K', 0.003578)
concrete.add_element('Ca', 0.014924)
concrete.add_element('Fe', 0.001605)
concrete.set_density('g/cm3', 2.3)

polyethylene = openmc.Material(name='Polyethylene, Borated')
polyethylene.add_element('H', 0.627756)
polyethylene.add_element('B', 0.046680)
polyethylene.add_element('C', 0.325564)
polyethylene.set_density('g/cm3', 1.0)

mats_conv = {
    mat.name: mat for mat in [steel, concrete, polyethylene]
}

openmc.Materials(mats_conv.values())

plotxs(mats_conv)

import pandas as pd

# Atomic fractions (ao).
compositions_ao = {
    'Steel, Stainless 316L': {
        'C': 0.001384,
        'Mn': 0.020165,
        'P': 0.000805,
        'S': 0.000518,
        'Si': 0.019722,
        'Cr': 0.181098,
        'Ni': 0.113247,
        'Mo': 0.014432,
        'Fe': 0.648629,
    },
    'Concrete, Ordinary (NIST)': {
        'H': 0.305322,
        'C': 0.002880,
        'O': 0.500411,
        'Na': 0.009212,
        'Mg': 0.000725,
        'Al': 0.010298,
        'Si': 0.151046,
        'K': 0.003578,
        'Ca': 0.014924,
        'Fe': 0.001605,
    },
    'Polyethylene, Borated': {
        'H': 0.627756,
        'B': 0.046680,
        'C': 0.325564,
    },
}

# Use column names matching the arguments of openmc.Material.add_element():
# https://docs.openmc.org/en/stable/pythonapi/generated/openmc.Material.html
# The next cell passes these as keyword arguments. For example,
# renaming 'percent' to 'fraction' would raise an error:
# TypeError: Material.add_element() got an unexpected keyword
# argument 'fraction'.
mats_df = pd.DataFrame(columns=['name', 'element', 'percent'])

for name, composition_ao in compositions_ao.items():
    for element, fraction in composition_ao.items():
        mats_df.loc[len(mats_df)] = [name, element, fraction]

mats_df

# Define the elemental composition.
mats_prog = {}
for name, dat in mats_df.groupby('name'):
    mats_prog[name] = openmc.Material(name=name)

    for _, row in dat.drop(columns='name').iterrows():
        mats_prog[name].add_element(**row.dropna().to_dict())

# Define the density.
densities = {
    'Steel, Stainless 316L': 8.0,
    'Concrete, Ordinary (NIST)': 2.3,
    'Polyethylene, Borated': 1.0,
}

for name, density in densities.items():
    mats_prog[name].set_density('g/cm3', density)
openmc.Materials(mats_prog.values())
# Three OpenMC materials successfully created.

plotxs(mats_prog)

import json
import subprocess

pnnl_url = (
    'https://raw.githubusercontent.com/openmc-dev/openmc/'
    'develop/openmc/data/material_libraries/pnnl_v2.json'
)
pnnl_file = DATA_ROOT / 'pnnl_v2.json'
download_pnnl = True  # Set False to reuse the downloaded file.

if download_pnnl:
    subprocess.run(
        ['curl', '-fLsS', '-o', str(pnnl_file), pnnl_url],
        check=True,
    )

if not pnnl_file.is_file():
    raise FileNotFoundError(
        f'{pnnl_file} is missing. Set download_pnnl = True.'
    )

with pnnl_file.open() as f:
    pnnl = json.load(f)

mats_pnnl = {}

for name in compositions_ao:
    data = pnnl['materials'][name]
    mat = openmc.Material(name=name)

    for element, fraction in data['elements'].items():
        mat.add_element(
            element,
            fraction,
            percent_type=pnnl['percent_type'],
        )

    mat.set_density(pnnl['density_units'], data['density'])
    mats_pnnl[name] = mat

openmc.Materials(mats_pnnl.values())

plotxs(mats_pnnl)

import neutronics_material_maker as nmm

available = nmm.AvailableMaterials()
list(available)

mats_nmm = {}

for name in compositions_ao:
    mat = nmm.Material.from_library(name=name).openmc_material
    mats_nmm[name] = mat

openmc.Materials(mats_nmm.values())

plotxs(mats_nmm)

mats = ['lithium-lead', 'FLiBe', 'Lithium', 'H2O']
for mat in mats:
    for kk, vv in available[mat].items():
        print('{:20s} {}'.format(kk, vv))
    print()

for mat in mats:
    density = available[mat]['density']
    material = nmm.Material.from_library(
        name=mat,
        temperature=294,
        pressure=1e5,
    ).openmc_material
    mass_density = material.get_mass_density()

    print(mat)
    print(f'  NMM density ({type(density).__name__}): {density}')
    print(f'  At 294 K, 1 bar: {mass_density:.6f} g/cm³')
    print()