Merge pull request #1523 from Mikolaj-A-Kowalski/cam_WP7_task4

Deterministic Cell Inventory Calculations

docs/source/usersguide/geometry.rst

@@ -185,7 +185,13 @@ the :class:`openmc.Cell` class::
 
 In this example, an instance of :class:`openmc.Material` is assigned to the
 :attr:`Cell.fill` attribute. One can also fill a cell with a :ref:`universe
-<usersguide_universes>` or :ref:`lattice <usersguide_lattices>`.
+<usersguide_universes>` or :ref:`lattice <usersguide_lattices>`. If you provide
+no fill to a cell or assign a value of `None`, it will be treated as a "void"
+cell with no material within. Particles are allowed to stream through the cell but
+will undergo no collisions::
+
+  # This cell will be filled with void on export to XML
+  gap = openmc.Cell(region=pellet_gap)
 
 The classes :class:`Halfspace`, :class:`Intersection`, :class:`Union`, and
 :class:`Complement` and all instances of :class:`openmc.Region` and can be
@@ -434,3 +440,30 @@ named ``dagmc.h5m``) when initializing a simulation. If a `geometry.xml
   <https://svalinn.github.io/DAGMC/usersguide/tools.html#make-watertight>`_. Future
   implementations of DAGMC geometry will support small volume overlaps and
   un-merged surfaces.
+
+-------------------------
+Calculating Atoms Content
+-------------------------
+
+If the total volume occupied by all instances of a cell in the geometry is known
+by the user, it is possible to assign this volume to a cell without performing a
+:ref:`stochastic volume <usersguide_volume>` calculation::
+
+  from uncertainties import ufloat
+
+  # Set known total volume in [cc]
+  cell = openmc.Cell()
+  cell.volume = 17.0
+
+  # Set volume if it is known with some uncertainty
+  cell.volume = ufloat(17.0, 0.1)
+
+Once a volume is set, and a cell is filled with a material or distributed
+materials, it is possible to use the :func:`~openmc.Cell.atoms` method to obtain
+a dictionary of nuclides and their total number of atoms in all instances
+of a cell (e.g. ``{'H1': 1.0e22, 'O16': 0.5e22, ...}``)::
+
+  cell = openmc.Cell(fill = u02)
+  cell.volume = 17.0
+
+  O16_atoms = cell.atoms['O16']

openmc/cell.py

@@ -4,6 +4,7 @@ from copy import deepcopy
 from math import cos, sin, pi
 from numbers import Real
 from xml.etree import ElementTree as ET
+from uncertainties import UFloat
 import sys
 import warnings
 
@@ -85,7 +86,12 @@ class Cell(IDManagerMixin):
     volume : float
         Volume of the cell in cm^3. This can either be set manually or
         calculated in a stochastic volume calculation and added via the
-        :meth:`Cell.add_volume_information` method.
+        :meth:`Cell.add_volume_information` method. For 'distribmat' cells
+        it is the total volume of all instances.
+    atoms : collections.OrderedDict
+        Mapping of nuclides to the total number of atoms for each nuclide
+        present in the cell, or in all of its instances for a 'distribmat'
+        fill. For example, {'U235': 1.0e22, 'U238': 5.0e22, ...}.
 
     """
 
@@ -134,6 +140,7 @@ class Cell(IDManagerMixin):
             string += '\t{0: <15}=\t{1}\n'.format('Temperature',
                                                   self.temperature)
         string += '{: <16}=\t{}\n'.format('\tTranslation', self.translation)
+        string += '{: <16}=\t{}\n'.format('\tVolume', self.volume)
 
         return string
 
@@ -182,6 +189,55 @@ class Cell(IDManagerMixin):
     def volume(self):
         return self._volume
 
+    @property
+    def atoms(self):
+        if self._atoms is None:
+            if self._volume is None:
+                msg = ('Cannot calculate atom content becouse no volume '
+                       'is set. Use Cell.volume to provide it or perform '
+                       'a stochastic volume calculation.')
+                raise ValueError(msg)
+
+            elif self.fill_type == 'void':
+                msg = ('Cell is filled with void. It contains no atoms. '
+                       'Material must be set to calculate atom content.')
+                raise ValueError(msg)
+
+            elif self.fill_type in ['lattice', 'universe']:
+                msg = ('Universe and Lattice cells can contain multiple '
+                       'materials in diffrent proportions. Atom content must '
+                       'be calculated with stochastic volume calculation.')
+                raise ValueError(msg)
+
+            elif self.fill_type == 'material':
+                # Get atomic densities
+                self._atoms = self._fill.get_nuclide_atom_densities()
+
+                # Convert to total number of atoms
+                for key, nuclide in self._atoms.items():
+                    atom = nuclide[1] * self._volume * 1.0e+24
+                    self._atoms[key] = atom
+
+            elif self.fill_type == 'distribmat':
+                # Assumes that volume is total volume of all instances
+                # Also assumes that all instances have the same volume
+                partial_volume = self.volume / len(self.fill)
+                self._atoms = OrderedDict()
+                for mat in self.fill:
+                    for key, nuclide in mat.get_nuclide_atom_densities().items():
+                        # To account for overlap of nuclides between distribmat
+                        # we need to append new atoms to any existing value
+                        # hence it is necessary to ask for default.
+                        atom = self._atoms.setdefault(key, 0)
+                        atom += nuclide[1] * partial_volume * 1.0e+24
+                        self._atoms[key] = atom
+
+            else:
+                msg = 'Unrecognized fill_type: {}'.format(self.fill_type)
+                raise ValueError(msg)
+
+        return self._atoms
+
     @property
     def paths(self):
         if self._paths is None:
@@ -223,12 +279,16 @@ class Cell(IDManagerMixin):
 
             elif not isinstance(fill, (openmc.Material, openmc.Lattice,
                                        openmc.Universe)):
-                msg = 'Unable to set Cell ID="{0}" to use a non-Material or ' \
-                      'Universe fill "{1}"'.format(self._id, fill)
+                msg = ('Unable to set Cell ID="{0}" to use a non-Material or '
+                       'Universe fill "{1}"'.format(self._id, fill))
                 raise ValueError(msg)
 
         self._fill = fill
 
+        # Info about atom content can now be invalid
+        # (since fill has just changed)
+        self._atoms = None
+
     @rotation.setter
     def rotation(self, rotation):
         cv.check_length('cell rotation', rotation, 3)
@@ -285,9 +345,15 @@ class Cell(IDManagerMixin):
     @volume.setter
     def volume(self, volume):
         if volume is not None:
-            cv.check_type('cell volume', volume, Real)
+            cv.check_type('cell volume', volume, (Real, UFloat))
+            cv.check_greater_than('cell volume', volume, 0.0, equality=True)
+
         self._volume = volume
 
+        # Info about atom content can now be invalid
+        # (sice volume has just changed)
+        self._atoms = None
+
     def add_volume_information(self, volume_calc):
         """Add volume information to a cell.
 

tests/regression_tests/distribmat/results_true.dat

@@ -7,3 +7,4 @@ Cell
 	Region         =	-9
 	Rotation       =	None
 	Translation    =	None
+	Volume         =	None

tests/regression_tests/multipole/results_true.dat

@@ -39,3 +39,4 @@ Cell
 	Rotation       =	None
 	Temperature    =	[500. 700.   0. 800.]
 	Translation    =	None
+	Volume         =	None

tests/unit_tests/test_cell.py

@@ -1,10 +1,13 @@
 import xml.etree. ElementTree as ET
 
 import numpy as np
+from uncertainties import ufloat
 import openmc
 import pytest
 
+
 from tests.unit_tests import assert_unbounded
+from openmc.data import atomic_mass, AVOGADRO
 
 
 def test_contains():
@@ -29,6 +32,14 @@ def test_repr(cell_with_lattice):
     c = openmc.Cell()
     repr(c)
 
+    # Empty cell with volume
+    c.volume = 3.0
+    repr(c)
+
+    # Empty cell with uncertain volume
+    c.volume = ufloat(3.0, 0.2)
+    repr(c)
+
 
 def test_bounding_box():
     zcyl = openmc.ZCylinder()
@@ -112,6 +123,114 @@ def test_get_nuclides(uo2):
     assert nucs == ['U235', 'O16']
 
 
+def test_volume_setting():
+    c = openmc.Cell()
+
+    # Test ordinary volume and uncertain volume
+    c.volume = 3
+    c.volume = ufloat(3, 0.7)
+
+    # Allow volume to be set to 0.0
+    c.volume = 0.0
+    c.volume = ufloat(0.0, 0.1)
+
+    # Test errors for negative volume
+    with pytest.raises(ValueError):
+        c.volume = -1.0
+    with pytest.raises(ValueError):
+        c.volume = ufloat(-0.05, 0.1)
+
+
+def test_atoms_material_cell(uo2, water):
+    """ Test if correct number of atoms is returned.
+    Also check if Cell.atoms still works after volume/material was changed
+    """
+    c = openmc.Cell(fill=uo2)
+    c.volume = 2.0
+    expected_nucs = ['U235', 'O16']
+
+    # Precalculate the expected number of atoms
+    M = (atomic_mass('U235') + 2 * atomic_mass('O16')) / 3
+    expected_atoms = [
+        1/3 * uo2.density/M * AVOGADRO * 2.0,  # U235
+        2/3 * uo2.density/M * AVOGADRO * 2.0   # O16
+    ]
+
+    tuples = c.atoms.items()
+    for nuc, atom_num, t in zip(expected_nucs, expected_atoms, tuples):
+        assert nuc == t[0]
+        assert atom_num == t[1]
+
+    # Change volume and check if OK
+    c.volume = 3.0
+    expected_atoms = [
+        1/3 * uo2.density/M * AVOGADRO * 3.0,  # U235
+        2/3 * uo2.density/M * AVOGADRO * 3.0   # O16
+    ]
+
+    tuples = c.atoms.items()
+    for nuc, atom_num, t in zip(expected_nucs, expected_atoms, tuples):
+        assert nuc == t[0]
+        assert atom_num == pytest.approx(t[1])
+
+    # Change material and check if OK
+    c.fill = water
+    expected_nucs = ['H1', 'O16']
+    M = (2 * atomic_mass('H1') + atomic_mass('O16')) / 3
+    expected_atoms = [
+        2/3 * water.density/M * AVOGADRO * 3.0,  # H1
+        1/3 * water.density/M * AVOGADRO * 3.0   # O16
+    ]
+
+    tuples = c.atoms.items()
+    for nuc, atom_num, t in zip(expected_nucs, expected_atoms, tuples):
+        assert nuc == t[0]
+        assert atom_num == pytest.approx(t[1])
+
+
+def test_atoms_distribmat_cell(uo2, water):
+    """ Test if correct number of atoms is returned for a cell with
+    'distribmat' fill
+    """
+    c = openmc.Cell(fill=[uo2, water])
+    c.volume = 6.0
+
+    # Calculate the expected number of atoms
+    expected_nucs = ['U235', 'O16', 'H1']
+    M_uo2 = (atomic_mass('U235') + 2 * atomic_mass('O16')) / 3
+    M_water = (2 * atomic_mass('H1') + atomic_mass('O16')) / 3
+    expected_atoms = [
+        1/3 * uo2.density/M_uo2 * AVOGADRO * 3.0,        # U235
+        (2/3 * uo2.density/M_uo2 * AVOGADRO * 3.0 +
+         1/3 * water.density/M_water * AVOGADRO * 3.0),  # O16
+        2/3 * water.density/M_water * AVOGADRO * 3.0     # H1
+    ]
+
+    tuples = c.atoms.items()
+    for nuc, atom_num, t in zip(expected_nucs, expected_atoms, tuples):
+        assert nuc == t[0]
+        assert atom_num == pytest.approx(t[1])
+
+
+def test_atoms_errors(cell_with_lattice):
+    cells, mats, univ, lattice = cell_with_lattice
+
+    # Material Cell with no volume
+    with pytest.raises(ValueError):
+        cells[1].atoms
+
+    # Cell with lattice
+    cells[2].volume = 3
+    with pytest.raises(ValueError):
+        cells[2].atoms
+
+    # Cell with volume but with void fill
+    cells[1].volume = 2
+    cells[1].fill = None
+    with pytest.raises(ValueError):
+        cells[1].atoms
+
+
 def test_nuclide_densities(uo2):
     c = openmc.Cell(fill=uo2)
     expected_nucs = ['U235', 'O16']
@@ -119,7 +238,7 @@ def test_nuclide_densities(uo2):
     tuples = list(c.get_nuclide_densities().values())
     for nuc, density, t in zip(expected_nucs, expected_density, tuples):
         assert nuc == t[0]
-        assert density == t[1]
+        assert density == pytest.approx(t[1])
 
     # Empty cell
     c = openmc.Cell()