From a06190ee407c14fd7a82cd0edf2fd89ce7e956d1 Mon Sep 17 00:00:00 2001
From: Sterling Harper <smharper@mit.edu>
Date: Sat, 10 Feb 2018 02:23:50 -0500
Subject: [PATCH] Add convenience function for boric acid Materials

---
 openmc/material.py                | 70 +++++++++++++++++++++++++++++++
 tests/unit_tests/test_material.py | 18 ++++++++
 2 files changed, 88 insertions(+)

diff --git a/openmc/material.py b/openmc/material.py
index e409d6536c..32b7d3a80d 100644
--- a/openmc/material.py
+++ b/openmc/material.py
@@ -973,3 +973,73 @@ class Materials(cv.CheckedList):
         # Write the XML Tree to the materials.xml file
         tree = ET.ElementTree(root_element)
         tree.write(path, xml_declaration=True, encoding='utf-8')
+
+
+def make_boric_acid(boron_ppm, temperature=293., pressure=0.1013, density=None,
+                    **kwargs):
+    """Return a Material with the composition of boric acid.
+
+    The water density can either be given directly, or it can be determined from
+    a temperature and pressure.
+
+    Parameters
+    ----------
+    boron_ppm : float
+        The weight fraction in parts-per-million of elemental boron in the acid.
+    temperature : float
+        Water temperature in [K] used to compute water density.
+    pressure : float
+        Water pressure in [MPa] used to compute water density.
+    density : float
+        Water density in [g / cm^3].  If specified, this value overrides the
+        temperature and pressure arguments.
+    **kwargs
+        All keyword arguments are passed to the created Material object.
+
+    Returns
+    -------
+    openmc.Material
+
+    """
+    # Set the density of water, either from an explicitly given density or from
+    # temperature and pressure.
+    if density is not None:
+        water_density = density
+    else:
+        water_density = openmc.data.water_density(temperature, pressure)
+
+    # Compute the density of the boric acid.
+    acid_density = water_density / (1 - boron_ppm * 1e-6)
+
+    # Compute the molar mass of pure water.
+    hydrogen = openmc.Element('H')
+    oxygen = openmc.Element('O')
+    M_H2O = 0.0
+    for iso_name, frac, junk in hydrogen.expand(2.0, 'ao'):
+        M_H2O += frac * openmc.data.atomic_mass(iso_name)
+    for iso_name, frac, junk in oxygen.expand(1.0, 'ao'):
+        M_H2O += frac * openmc.data.atomic_mass(iso_name)
+
+    # Compute the molar mass of boron.
+    boron = openmc.Element('B')
+    M_B = 0.0
+    for iso_name, frac, junk in boron.expand(1.0, 'ao'):
+        M_B += frac * openmc.data.atomic_mass(iso_name)
+
+    # Compute the number fractions of each element.
+    frac_H2O = (1 - boron_ppm * 1e-6) / M_H2O
+    frac_H = 2 * frac_H2O
+    frac_O = frac_H2O
+    frac_B = boron_ppm * 1e-6 / M_B
+
+    # Build the material.
+    if density is None:
+        out = openmc.Material(temperature=temperature, **kwargs)
+    else:
+        out = openmc.Material(**kwargs)
+    out.add_element('H', frac_H, 'ao')
+    out.add_element('O', frac_O, 'ao')
+    out.add_element('B', frac_B, 'ao')
+    out.set_density('g/cc', acid_density)
+    out.add_s_alpha_beta('c_H_in_H2O')
+    return out
diff --git a/tests/unit_tests/test_material.py b/tests/unit_tests/test_material.py
index e92a2ac088..38495cb1ab 100644
--- a/tests/unit_tests/test_material.py
+++ b/tests/unit_tests/test_material.py
@@ -139,3 +139,21 @@ def test_materials(run_in_tmpdir):
     mats.cross_sections = '/some/fake/cross_sections.xml'
     mats.multipole_library = '/some/awesome/mp_lib/'
     mats.export_to_xml()
+
+
+def test_boric_acid():
+    # Test against reference values from the BEAVRS benchmark.
+    m = openmc.make_boric_acid(975, 566.5, 15.51, material_id=50)
+    assert m.density == pytest.approx(0.7405, 1e-3)
+    assert m.temperature == pytest.approx(566.5)
+    assert m._sab[0][0] == 'c_H_in_H2O'
+    ref_dens = {'B10':8.0023e-06, 'B11':3.2210e-05, 'H1':4.9458e-02,
+                'O16':2.4672e-02}
+    nuc_dens = m.get_nuclide_atom_densities()
+    for nuclide in ref_dens:
+        assert nuc_dens[nuclide][1] == pytest.approx(ref_dens[nuclide], 1e-2)
+    assert m.id == 50
+
+    # Make sure the density override works
+    m = openmc.make_boric_acid(975, 566.5, 15.51, 0.9)
+    assert m.density == pytest.approx(0.9, 1e-3)