Repository: openmc-dev/openmc

Issue:
Extra error checking on mesh creation
## Description
Currently it is possible to build cylindricalmesh objects with phi_grid values that don't make sense. For example

```python
cylindrical_mesh = openmc.CylindricalMesh(
    r_grid=[0, 1, 2],
    phi_grid=[2],
    z_grid=[0, 10, 20],
)
>>> CylindricalMesh
        ID             =        9
        Name           =
        Dimensions     =        3
        Origin         =        [0. 0. 0.]
        N R pnts:      =        3
        R Min:         =        0.0
        R Max:         =        1.0
        N Phi pnts:    =        1
        Phi Min:       =        2.0   <----- min is the same as max
        Phi Max:       =        2.0   <----- min is the same as max
        N Z pnts:      =        4
        Z Min:         =        0.0
        Z Max:         =        20.0
```

Perhaps we could add some more checking on the setters to check that the min value and max value for phi are not the same and raise a ValueError if they are the same. I guess this is the same for r_grid and z_grid values and even for other mesh types.

## Alternatives
We could leave the code as is, The error gets reported to the user when they run the openmc executable. I think reporting the error when running the simulation is less desirable and Ideally the user would find out earlier (when creating the mesh)

## Compatibility
no API changes needed


Discussion:
hsameer481: @shimwell We have made the following changes to the cylindrical mesh based on the https://github.com/openmc-dev/openmc/issues/2933 issue. As first-time (Not just OpenMC but also GitHub) contributors, we wanted to ask if you could review and recommend any changes. 
shimwell: Wow this is some nice input checking, I hadn't though of checking the numbers increase.

It is certainly an improvement on what we have at the moment.

Would the increasing check catch an input like ```r_grid=[1,2,4,3]```
hsameer481: @shimwell Thank you for getting back to us! We have committed new changes that account for input checks similar to [1, 2, 4, 3] as well as adding a unit test for this. Does this logic make sense?
shimwell: This all looks good to me, there [are lots of ways of checking the list is ascending](https://www.geeksforgeeks.org/python-check-if-list-is-strictly-increasing/) so I guess we want to make sure we have a fast one.

Also I think this is useful enough in multiple places that we could consider making a check_ascending function in [this file ](https://github.com/openmc-dev/openmc/blob/develop/openmc/checkvalue.py). and making use of it in a few places throughout the code.

As a user who makes tons of mistakes I'm always keen to see extra input checking and am in favour of the PR and also like the code and tests you have done.

I'm going to defer to @paulromano for the final judgment on this PR to check we have the balance of extra compute with benefit to user correct.
shimwell: @hsameer481 what do you think about the code in #2973. The code added to checkvalues.py looks like we can make use of it on your PR.

What do you think should we merge in #2973 first, then sync your fork and make use of those functions on this PR
hsameer481: @shimwell that sounds good to us! Since this is the first time we're doing this, is there a specific process for implementing a new merged PR? 
shimwell: > @shimwell that sounds good to us! Since this is the first time we're doing this, is there a specific process for implementing a new merged PR?

Yep, once the PR has been merged (not yet)

1. go to you own fork
2. selected develop branch on the dropdown of branches
3. click "sync fork"
4. go to your local computer terminal and cd into the repository folder
5. git checkout develop
6. git pull
7. git checkout new_cylindrical_mesh
8. git merge develop
9. then you should have access to the new code which was pushed to openmc-dev/openmc:develop on your feature branch. So you could make use of the new check code.

But we must first wait to see the other PR gets merged
hsameer481: @shimwell Great! We really appreciate the guidance and will make sure to follow those steps once it's merged!

Patch:
diff --git a/openmc/lattice.py b/openmc/lattice.py
index f7f9da8ea..518068560 100644
--- a/openmc/lattice.py
+++ b/openmc/lattice.py
@@ -884,6 +884,10 @@ class RectLattice(Lattice):
         dimension = ET.SubElement(lattice_subelement, "dimension")
         dimension.text = ' '.join(map(str, self.shape))
 
+        # Make sure lower_left has been specified
+        if self.lower_left is None:
+            raise ValueError(f"Lattice {self.id} does not have lower_left specified.")
+
         # Export Lattice lower left
         lower_left = ET.SubElement(lattice_subelement, "lower_left")
         lower_left.text = ' '.join(map(str, self._lower_left))
diff --git a/openmc/mesh.py b/openmc/mesh.py
index 48263c205..215c91923 100644
--- a/openmc/mesh.py
+++ b/openmc/mesh.py
@@ -152,6 +152,7 @@ class MeshBase(IDManagerMixin, ABC):
             model: openmc.Model,
             n_samples: int = 10_000,
             prn_seed: Optional[int] = None,
+            include_void: bool = True,
             **kwargs
     ) -> List[openmc.Material]:
         """Generate homogenized materials over each element in a mesh.
@@ -168,6 +169,8 @@ class MeshBase(IDManagerMixin, ABC):
         prn_seed : int, optional
             Pseudorandom number generator (PRNG) seed; if None, one will be
             generated randomly.
+        include_void : bool, optional
+            Whether homogenization should include voids.
         **kwargs
             Keyword-arguments passed to :func:`openmc.lib.init`.
 
@@ -222,6 +225,10 @@ class MeshBase(IDManagerMixin, ABC):
                 material_ids.pop(index_void)
                 volumes.pop(index_void)
 
+            # If void should be excluded, adjust total volume
+            if not include_void:
+                total_volume = sum(volumes)
+
             # Compute volume fractions
             volume_fracs = np.array(volumes) / total_volume
 
@@ -1381,6 +1388,8 @@ class CylindricalMesh(StructuredMesh):
     @r_grid.setter
     def r_grid(self, grid):
         cv.check_type('mesh r_grid', grid, Iterable, Real)
+        cv.check_length('mesh r_grid', grid, 2)
+        cv.check_increasing('mesh r_grid', grid)
         self._r_grid = np.asarray(grid, dtype=float)
 
     @property
@@ -1390,7 +1399,12 @@ class CylindricalMesh(StructuredMesh):
     @phi_grid.setter
     def phi_grid(self, grid):
         cv.check_type('mesh phi_grid', grid, Iterable, Real)
-        self._phi_grid = np.asarray(grid, dtype=float)
+        cv.check_length('mesh phi_grid', grid, 2)
+        cv.check_increasing('mesh phi_grid', grid)
+        grid = np.asarray(grid, dtype=float)
+        if np.any((grid < 0.0) | (grid > 2*pi)):
+            raise ValueError("phi_grid values must be in [0, 2π].")
+        self._phi_grid = grid
 
     @property
     def z_grid(self):
@@ -1399,6 +1413,8 @@ class CylindricalMesh(StructuredMesh):
     @z_grid.setter
     def z_grid(self, grid):
         cv.check_type('mesh z_grid', grid, Iterable, Real)
+        cv.check_length('mesh z_grid', grid, 2)
+        cv.check_increasing('mesh z_grid', grid)
         self._z_grid = np.asarray(grid, dtype=float)
 
     @property
@@ -1833,7 +1849,10 @@ class SphericalMesh(StructuredMesh):
         cv.check_type('mesh theta_grid', grid, Iterable, Real)
         cv.check_length('mesh theta_grid', grid, 2)
         cv.check_increasing('mesh theta_grid', grid)
-        self._theta_grid = np.asarray(grid, dtype=float)
+        grid = np.asarray(grid, dtype=float)
+        if np.any((grid < 0.0) | (grid > pi)):
+            raise ValueError("theta_grid values must be in [0, π].")
+        self._theta_grid = grid
 
     @property
     def phi_grid(self):
@@ -1844,7 +1863,10 @@ class SphericalMesh(StructuredMesh):
         cv.check_type('mesh phi_grid', grid, Iterable, Real)
         cv.check_length('mesh phi_grid', grid, 2)
         cv.check_increasing('mesh phi_grid', grid)
-        self._phi_grid = np.asarray(grid, dtype=float)
+        grid = np.asarray(grid, dtype=float)
+        if np.any((grid < 0.0) | (grid > 2*pi)):
+            raise ValueError("phi_grid values must be in [0, 2π].")
+        self._phi_grid = grid
 
     @property
     def _grids(self):
diff --git a/openmc/plotter.py b/openmc/plotter.py
index 97ca5f929..aa3825f5e 100644
--- a/openmc/plotter.py
+++ b/openmc/plotter.py
@@ -59,10 +59,15 @@ def _get_legend_label(this, type):
     """Gets a label for the element or nuclide or material and reaction plotted"""
     if isinstance(this, str):
         if type in openmc.data.DADZ:
-            z, a, m = openmc.data.zam(this)
-            da, dz = openmc.data.DADZ[type]
-            gnds_name = openmc.data.gnds_name(z + dz, a + da, m)
-            return f'{this} {type} {gnds_name}'
+            if this in ELEMENT_NAMES:
+                return f'{this} {type}'
+            else:  # this is a nuclide so the legend can contain more information
+                z, a, m = openmc.data.zam(this)
+                da, dz = openmc.data.DADZ[type]
+                gnds_name = openmc.data.gnds_name(z + dz, a + da, m)
+                # makes a string with nuclide reaction and new nuclide
+                # For example "Be9 (n,2n) Be8"
+                return f'{this} {type} {gnds_name}'
         return f'{this} {type}'
     elif this.name == '':
         return f'Material {this.id} {type}'


Test patch:
diff --git a/tests/unit_tests/mesh_to_vtk/test_vtk_dims.py b/tests/unit_tests/mesh_to_vtk/test_vtk_dims.py
index baa544984..0b6acf5aa 100644
--- a/tests/unit_tests/mesh_to_vtk/test_vtk_dims.py
+++ b/tests/unit_tests/mesh_to_vtk/test_vtk_dims.py
@@ -153,7 +153,7 @@ def test_write_data_to_vtk_round_trip(run_in_tmpdir):
 
     smesh = openmc.SphericalMesh(
         r_grid=(0.0, 1.0, 2.0),
-        theta_grid=(0.0, 2.0, 4.0, 5.0),
+        theta_grid=(0.0, 0.5, 1.0, 2.0),
         phi_grid=(0.0, 3.0, 6.0),
     )
     rmesh = openmc.RegularMesh()
diff --git a/tests/unit_tests/test_mesh.py b/tests/unit_tests/test_mesh.py
index 5bacd5fc6..ed08be816 100644
--- a/tests/unit_tests/test_mesh.py
+++ b/tests/unit_tests/test_mesh.py
@@ -189,6 +189,42 @@ def test_CylindricalMesh_initiation():
         openmc.SphericalMesh(('🧇', '🥞'))
 
 
+def test_invalid_cylindrical_mesh_errors():
+    # Test invalid r_grid values
+    with pytest.raises(ValueError):
+        openmc.CylindricalMesh(r_grid=[5, 1], phi_grid=[0, pi], z_grid=[0, 10])
+
+    with pytest.raises(ValueError):
+        openmc.CylindricalMesh(r_grid=[1, 2, 4, 3], phi_grid=[0, pi], z_grid=[0, 10])
+
+    with pytest.raises(ValueError):
+        openmc.CylindricalMesh(r_grid=[1], phi_grid=[0, pi], z_grid=[0, 10])
+
+    # Test invalid phi_grid values
+    with pytest.raises(ValueError):
+        openmc.CylindricalMesh(r_grid=[0, 1, 2], phi_grid=[-1, 3], z_grid=[0, 10])
+
+    with pytest.raises(ValueError):
+        openmc.CylindricalMesh(
+            r_grid=[0, 1, 2],
+            phi_grid=[0, 2*pi + 0.1],
+            z_grid=[0, 10]
+        )
+
+    with pytest.raises(ValueError):
+        openmc.CylindricalMesh(r_grid=[0, 1, 2], phi_grid=[pi], z_grid=[0, 10])
+
+    # Test invalid z_grid values
+    with pytest.raises(ValueError):
+        openmc.CylindricalMesh(r_grid=[0, 1, 2], phi_grid=[0, pi], z_grid=[5])
+
+    with pytest.raises(ValueError):
+        openmc.CylindricalMesh(r_grid=[0, 1, 2], phi_grid=[0, pi], z_grid=[5, 1])
+
+    with pytest.raises(ValueError):
+        openmc.CylindricalMesh(r_grid=[1, 2, 4, 3], phi_grid=[0, pi], z_grid=[0, 10, 5])
+
+
 def test_centroids():
     # regular mesh
     mesh = openmc.RegularMesh()
@@ -386,3 +422,9 @@ def test_mesh_get_homogenized_materials():
 
     # Mesh element that overlaps void should have half density
     assert m4.get_mass_density('H1') == pytest.approx(0.5, rel=1e-2)
+
+    # If not including void, density of homogenized material should be same as
+    # original material
+    m5, = mesh_void.get_homogenized_materials(
+        model, n_samples=1000, include_void=False)
+    assert m5.get_mass_density('H1') == pytest.approx(1.0)
diff --git a/tests/unit_tests/test_plotter.py b/tests/unit_tests/test_plotter.py
index 2c195c5e1..0220cec3e 100644
--- a/tests/unit_tests/test_plotter.py
+++ b/tests/unit_tests/test_plotter.py
@@ -75,7 +75,7 @@ def test_calculate_cexs_with_materials(test_mat):
 @pytest.mark.parametrize("this", ["Be", "Be9"])
 def test_plot_xs(this):
     from matplotlib.figure import Figure
-    assert isinstance(openmc.plot_xs({this: ['total', 'elastic']}), Figure)
+    assert isinstance(openmc.plot_xs({this: ['total', 'elastic', 16, '(n,2n)']}), Figure)
 
 
 def test_plot_xs_mat(test_mat):