Merge branch 'develop' into centre_for_cylinder_spherical_meshes

tests/unit_tests/test_deplete_resultslist.py

@@ -22,7 +22,7 @@ def test_get_atoms(res):
 
     t_ref = np.array([0.0, 1296000.0, 2592000.0, 3888000.0])
     n_ref = np.array(
-        [6.67473282e+08, 3.72442707e+14, 3.61129692e+14, 4.01920099e+14])
+        [6.67473282e+08, 3.88942731e+14, 3.73091215e+14, 3.26987387e+14])
 
     np.testing.assert_allclose(t, t_ref)
     np.testing.assert_allclose(n, n_ref)
@@ -48,8 +48,8 @@ def test_get_reaction_rate(res):
     t, r = res.get_reaction_rate("1", "Xe135", "(n,gamma)")
 
     t_ref = [0.0, 1296000.0, 2592000.0, 3888000.0]
-    n_ref = [6.67473282e+08, 3.72442707e+14, 3.61129692e+14, 4.01920099e+14]
-    xs_ref = [5.10301159e-05, 3.19379638e-05, 4.50543806e-05, 4.71004301e-05]
+    n_ref = [6.67473282e+08, 3.88942731e+14, 3.73091215e+14, 3.26987387e+14]
+    xs_ref = [2.53336104e-05, 4.21747011e-05, 3.48616127e-05, 3.61775563e-05]
 
     np.testing.assert_allclose(t, t_ref)
     np.testing.assert_allclose(r, np.array(n_ref) * xs_ref)
@@ -61,8 +61,8 @@ def test_get_keff(res):
     t_min, k = res.get_keff(time_units='min')
 
     t_ref = [0.0, 1296000.0, 2592000.0, 3888000.0]
-    k_ref = [1.21409662, 1.16518654, 1.25357797, 1.22611968]
-    u_ref = [0.0278795195, 0.0233141097, 0.0167899218, 0.0246734716]
+    k_ref = [1.1596402556, 1.1914183335, 1.2292570871, 1.1797030302]
+    u_ref = [0.0270680649, 0.0219163444, 0.024268508 , 0.0221401194]
 
     np.testing.assert_allclose(t, t_ref)
     np.testing.assert_allclose(t_min * 60, t_ref)

tests/unit_tests/test_filter_mesh.py

@@ -1,7 +1,10 @@
+import math
+
 import numpy as np
-import openmc
 from uncertainties import unumpy
 
+import openmc
+
 
 def test_spherical_mesh_estimators(run_in_tmpdir):
     """Test that collision/tracklength estimators agree for SphericalMesh"""
@@ -62,7 +65,8 @@ def test_cylindrical_mesh_estimators(run_in_tmpdir):
     mat.add_nuclide('U235', 1.0)
     mat.set_density('g/cm3', 10.0)
 
-    cyl = openmc.model.RightCircularCylinder((0., 0., -5.), 10., 10.0, boundary_type='vacuum')
+    cyl = openmc.model.RightCircularCylinder((0., 0., -5.), 10., 10.0,
+                                             boundary_type='vacuum')
     cell = openmc.Cell(fill=mat, region=-cyl)
     model = openmc.Model()
     model.geometry = openmc.Geometry([cell])
@@ -107,3 +111,43 @@ def test_cylindrical_mesh_estimators(run_in_tmpdir):
     diff = unumpy.nominal_values(delta)
     std_dev = unumpy.std_devs(delta)
     assert np.all(diff < 3*std_dev)
+
+def test_get_reshaped_data(run_in_tmpdir):
+    """Test that expanding MeshFilter dimensions works as expected"""
+
+    mat = openmc.Material()
+    mat.add_nuclide('U235', 1.0)
+    mat.set_density('g/cm3', 10.0)
+
+    sphere = openmc.Sphere(r=10.0, boundary_type='vacuum')
+    cell = openmc.Cell(fill=mat, region=-sphere)
+    model = openmc.Model()
+    model.geometry = openmc.Geometry([cell])
+    model.settings.particles = 1_000
+    model.settings.inactive = 10
+    model.settings.batches = 20
+
+    sph_mesh = openmc.SphericalMesh()
+    sph_mesh.r_grid = np.linspace(0.0, 5.0**3, 20)**(1/3)
+    sph_mesh.theta_grid = np.linspace(0, math.pi, 4)
+    sph_mesh.phi_grid = np.linspace(0, 2*math.pi, 3)
+    tally1 = openmc.Tally()
+    efilter = openmc.EnergyFilter([0, 1e5, 1e8])
+    meshfilter = openmc.MeshFilter(sph_mesh)
+    assert meshfilter.shape == (19, 3, 2)
+    tally1.filters = [efilter, meshfilter]
+    tally1.scores = ['flux']
+
+    model.tallies = openmc.Tallies([tally1])
+
+    # Run OpenMC
+    sp_filename = model.run()
+
+    # Get flux tally as reshaped data
+    with openmc.StatePoint(sp_filename) as sp:
+        t1 = sp.tallies[tally1.id]
+        data1 = t1.get_reshaped_data()
+        data2 = t1.get_reshaped_data(expand_dims=True)
+
+    assert data1.shape == (2, 19*3*2, 1, 1)
+    assert data2.shape == (2, 19, 3, 2, 1, 1)

tests/unit_tests/test_filters.py

@@ -269,3 +269,23 @@ def test_energyfunc():
     np.testing.assert_allclose(f.energy, new_f.energy)
     np.testing.assert_allclose(f.y, new_f.y)
     assert f.interpolation == new_f.interpolation
+
+
+def test_tabular_from_energyfilter():
+    efilter = openmc.EnergyFilter([0.0, 10.0, 20.0, 25.0])
+    tab = efilter.get_tabular(values=[5, 10, 10])
+
+    assert tab.x.tolist() == [0.0, 10.0, 20.0, 25.0]
+
+    # combination of different values passed into get_tabular and different
+    # width energy bins results in a doubling value for each p value
+    assert tab.p.tolist() == [0.02, 0.04, 0.08, 0.0]
+
+    # distribution should integrate to unity
+    assert tab.integral() == approx(1.0)
+
+    # 'histogram' is the default
+    assert tab.interpolation == 'histogram'
+
+    tab = efilter.get_tabular(values=np.array([10, 10, 5]), interpolation='linear-linear')
+    assert tab.interpolation == 'linear-linear'

tests/unit_tests/test_mesh_tets.e

@@ -0,0 +1 @@
+../regression_tests/unstructured_mesh/test_mesh_tets.e

tests/unit_tests/test_plotter.py

@@ -1,21 +1,27 @@
-import openmc
 import numpy as np
+import openmc
+import pytest
+from matplotlib.figure import Figure
 
 
-def test_calculate_cexs_elem_mat_sab():
-    """Checks that sab cross sections are included in the
-    _calculate_cexs_elem_mat method and have the correct shape"""
-
+@pytest.fixture(scope="module")
+def test_mat():
     mat_1 = openmc.Material()
     mat_1.add_element("H", 4.0, "ao")
     mat_1.add_element("O", 4.0, "ao")
     mat_1.add_element("C", 4.0, "ao")
+    return mat_1
 
-    mat_1.add_s_alpha_beta("c_C6H6")
-    mat_1.set_density("g/cm3", 0.865)
+
+def test_calculate_cexs_elem_mat_sab(test_mat):
+    """Checks that sab cross sections are included in the
+    _calculate_cexs_elem_mat method and have the correct shape"""
+
+    test_mat.add_s_alpha_beta("c_C6H6")
+    test_mat.set_density("g/cm3", 0.865)
 
     energy_grid, data = openmc.plotter._calculate_cexs_elem_mat(
-        mat_1,
+        test_mat,
         ["inelastic"],
         sab_name="c_C6H6",
     )
@@ -25,3 +31,57 @@ def test_calculate_cexs_elem_mat_sab():
     assert len(energy_grid) > 1
     assert len(data) == 1
     assert len(data[0]) == len(energy_grid)
+
+
+@pytest.mark.parametrize("this,data_type", [("Li", "element"), ("Li6", "nuclide")])
+def test_calculate_cexs_with_element(this, data_type):
+
+    # single type (reaction)
+    energy_grid, data = openmc.plotter.calculate_cexs(
+        this=this, data_type=data_type, types=[205]
+    )
+
+    assert isinstance(energy_grid, np.ndarray)
+    assert isinstance(data, np.ndarray)
+    assert len(energy_grid) > 1
+    assert len(data) == 1
+    assert len(data[0]) == len(energy_grid)
+
+    # two types (reaction)
+    energy_grid, data = openmc.plotter.calculate_cexs(
+        this=this, data_type=data_type, types=[2, "elastic"]
+    )
+
+    assert isinstance(energy_grid, np.ndarray)
+    assert isinstance(data, np.ndarray)
+    assert len(energy_grid) > 1
+    assert len(data) == 2
+    assert len(data[0]) == len(energy_grid)
+    assert len(data[0]) == len(energy_grid)
+    # reactions are both the same MT number 2 is elastic
+    assert np.array_equal(data[0], data[1])
+
+
+def test_calculate_cexs_with_materials(test_mat):
+    energy_grid, data = openmc.plotter.calculate_cexs(
+        this=test_mat, types=[205], data_type="material"
+    )
+
+    assert isinstance(energy_grid, np.ndarray)
+    assert isinstance(data, np.ndarray)
+    assert len(energy_grid) > 1
+    assert len(data) == 1
+    assert len(data[0]) == len(energy_grid)
+
+
+@pytest.mark.parametrize(("this,data_type"), [("Be", "element"), ("Be9", "nuclide")])
+def test_plot_xs(this, data_type):
+    assert isinstance(
+        openmc.plotter.plot_xs(this, data_type=data_type, types=["total"]), Figure
+    )
+
+
+def test_plot_xs_mat(test_mat):
+    assert isinstance(
+        openmc.plotter.plot_xs(test_mat, data_type="material", types=["total"]), Figure
+    )

tests/unit_tests/test_source_mesh.py

@@ -0,0 +1,203 @@
+from itertools import product
+from pathlib import Path
+from subprocess import call
+
+import pytest
+import numpy as np
+import openmc
+import openmc.lib
+
+from tests import cdtemp
+from tests.regression_tests import config
+
+
+TETS_PER_VOXEL = 12
+
+# This test uses a geometry file with cells that match a regular mesh. Each cell
+# in the geometry corresponds to 12 tetrahedra in the unstructured mesh file.
+@pytest.fixture
+def model():
+    openmc.reset_auto_ids()
+
+    ### Materials ###
+    materials = openmc.Materials()
+
+    water_mat = openmc.Material(name="water")
+    water_mat.add_nuclide("H1", 2.0)
+    water_mat.add_nuclide("O16", 1.0)
+    water_mat.set_density("atom/b-cm", 0.07416)
+    materials.append(water_mat)
+
+    ### Geometry ###
+    # This test uses a geometry file that resembles a regular mesh.
+    # 12 tets are used to match each voxel in the geometry.
+
+    # create a regular mesh that matches the superimposed mesh
+    regular_mesh = openmc.RegularMesh(mesh_id=10)
+    regular_mesh.lower_left = (-10, -10, -10)
+    regular_mesh.dimension = (10, 10, 10)
+    regular_mesh.width = (2, 2, 2)
+
+    root_cell, _ = regular_mesh.build_cells(bc=['vacuum']*6)
+
+    geometry = openmc.Geometry(root=[root_cell])
+
+    ### Settings ###
+    settings = openmc.Settings()
+    settings.run_mode = 'fixed source'
+    settings.particles = 100
+    settings.batches = 2
+
+    return openmc.model.Model(geometry=geometry,
+                              materials=materials,
+                              settings=settings)
+
+### Setup test cases ###
+param_values = (['libmesh', 'moab'], # mesh libraries
+                ['uniform', 'manual']) # Element weighting schemes
+
+test_cases = []
+for i, (lib, schemes) in enumerate(product(*param_values)):
+    test_cases.append({'library' : lib,
+                       'source_strengths' : schemes})
+
+def ids(params):
+    """Test naming function for clarity"""
+    return f"{params['library']}-{params['source_strengths']}"
+
+@pytest.mark.parametrize("test_cases", test_cases, ids=ids)
+def test_unstructured_mesh_sampling(model, request, test_cases):
+    # skip the test if the library is not enabled
+    if test_cases['library'] == 'moab' and not openmc.lib._dagmc_enabled():
+        pytest.skip("DAGMC (and MOAB) mesh not enabled in this build.")
+
+    if test_cases['library'] == 'libmesh' and not openmc.lib._libmesh_enabled():
+        pytest.skip("LibMesh is not enabled in this build.")
+
+    # setup mesh source ###
+    mesh_filename = Path(request.fspath).parent / "test_mesh_tets.e"
+    uscd_mesh = openmc.UnstructuredMesh(mesh_filename, test_cases['library'])
+
+    # subtract one to account for root cell produced by RegularMesh.build_cells
+    n_cells = len(model.geometry.get_all_cells()) - 1
+
+    # set source weights according to test case
+    if test_cases['source_strengths'] == 'uniform':
+        vol_norm = True
+        strengths = None
+    elif test_cases['source_strengths'] == 'manual':
+        vol_norm = False
+        # assign random weights
+        strengths = np.random.rand(n_cells*TETS_PER_VOXEL)
+
+    # create the spatial distribution based on the mesh
+    space = openmc.stats.MeshSpatial(uscd_mesh, strengths, vol_norm)
+
+    energy = openmc.stats.Discrete(x=[15.e+06], p=[1.0])
+    source = openmc.Source(space=space, energy=energy)
+    model.settings.source = source
+
+    with cdtemp([mesh_filename]):
+        model.export_to_xml()
+
+        n_measurements = 100
+        n_samples = 1000
+
+        cell_counts = np.zeros((n_cells, n_measurements))
+
+        # This model contains 1000 geometry cells. Each cell is a hex
+        # corresponding to 12 of the tets. This test runs 1000 samples. This
+        #  results in the following average for each cell
+        openmc.lib.init([])
+
+        # perform many sets of samples and track counts for each cell
+        for m in range(n_measurements):
+            sites = openmc.lib.sample_external_source(n_samples)
+            cells = [openmc.lib.find_cell(s.r) for s in sites]
+
+            for c in cells:
+                # subtract one from index to account for root cell
+                cell_counts[c[0]._index - 1, m] += 1
+
+        # make sure particle transport is successful
+        openmc.lib.run()
+        openmc.lib.finalize()
+
+        # normalize cell counts to get sampling frequency per particle
+        cell_counts /= n_samples
+
+        # get the mean and std. dev. of the cell counts
+        mean = cell_counts.mean(axis=1)
+        std_dev = cell_counts.std(axis=1)
+
+        if test_cases['source_strengths'] == 'uniform':
+            exp_vals = np.ones(n_cells) / n_cells
+        else:
+            # sum up the source strengths for each tet, these are the expected true mean
+            # of the sampling frequency for that cell
+            exp_vals = strengths.reshape(-1, 12).sum(axis=1) / sum(strengths)
+
+        diff = np.abs(mean - exp_vals)
+        assert((diff < 2*std_dev).sum() / diff.size >= 0.95)
+        assert((diff < 6*std_dev).sum() / diff.size >= 0.997)
+
+
+def test_strengths_size_failure(request, model):
+    # setup mesh source ###
+    mesh_filename = Path(request.fspath).parent / "test_mesh_tets.e"
+    uscd_mesh = openmc.UnstructuredMesh(mesh_filename, 'libmesh')
+
+    # intentionally incorrectly sized to trigger an error
+    n_cells = len(model.geometry.get_all_cells())
+    strengths = np.random.rand(n_cells*TETS_PER_VOXEL)
+
+    # create the spatial distribution based on the mesh
+    space = openmc.stats.MeshSpatial(uscd_mesh, strengths)
+
+    energy = openmc.stats.Discrete(x=[15.e+06], p=[1.0])
+    source = openmc.Source(space=space, energy=energy)
+    model.settings.source = source
+
+    # skip the test if unstructured mesh is not available
+    if not openmc.lib._libmesh_enabled():
+        if openmc.lib._dagmc_enabled():
+            source.space.mesh.library = 'moab'
+        else:
+            pytest.skip("Unstructured mesh support unavailable.")
+
+    # make sure that an incorrrectly sized strengths array causes a failure
+    source.space.strengths = source.space.strengths[:-1]
+
+    mesh_filename = Path(request.fspath).parent / source.space.mesh.filename
+
+    with pytest.raises(RuntimeError, match=r'strengths array'), cdtemp([mesh_filename]):
+        model.export_to_xml()
+        openmc.run()
+
+def test_roundtrip(run_in_tmpdir, model, request):
+    if not openmc.lib._libmesh_enabled() and not openmc.lib._dagmc_enabled():
+        pytest.skip("Unstructured mesh is not enabled in this build.")
+
+    mesh_filename = Path(request.fspath).parent / 'test_mesh_tets.e'
+    ucd_mesh = openmc.UnstructuredMesh(mesh_filename, library='libmesh')
+
+    if not openmc.lib._libmesh_enabled():
+        ucd_mesh.library = 'moab'
+
+    n_cells = len(model.geometry.get_all_cells())
+
+    space_out = openmc.MeshSpatial(ucd_mesh)
+    space_out.strengths = np.random.rand(n_cells*TETS_PER_VOXEL)
+    model.settings.source = openmc.Source(space=space_out)
+
+    # write out the model
+    model.export_to_xml()
+
+    model_in = openmc.Model.from_xml()
+
+    space_in = model_in.settings.source[0].space
+
+    np.testing.assert_equal(space_out.strengths, space_in.strengths)
+
+    assert space_in.mesh.id == space_out.mesh.id
+    assert space_in.volume_normalized == space_out.volume_normalized

tests/unit_tests/test_stats.py

@@ -166,7 +166,7 @@ def test_watt():
 
 def test_tabular():
     x = np.array([0.0, 5.0, 7.0])
-    p = np.array([0.1, 0.2, 0.05])
+    p = np.array([10.0, 20.0, 5.0])
     d = openmc.stats.Tabular(x, p, 'linear-linear')
     elem = d.to_xml_element('distribution')
 
@@ -178,19 +178,22 @@ def test_tabular():
 
     # test linear-linear sampling
     d = openmc.stats.Tabular(x, p)
-
     n_samples = 100_000
     samples = d.sample(n_samples)
     assert_sample_mean(samples, d.mean())
 
-    # test histogram sampling
-    d = openmc.stats.Tabular(x, p, interpolation='histogram')
+    # test linear-linear normalization
     d.normalize()
     assert d.integral() == pytest.approx(1.0)
 
+    # test histogram sampling
+    d = openmc.stats.Tabular(x, p, interpolation='histogram')
     samples = d.sample(n_samples)
     assert_sample_mean(samples, d.mean())
 
+    d.normalize()
+    assert d.integral() == pytest.approx(1.0)
+
 
 def test_legendre():
     # Pu239 elastic scattering at 100 keV

tests/unit_tests/test_universe.py

@@ -93,10 +93,12 @@ def test_get_all_universes():
     u2 = openmc.Universe(cells=[c2])
     c3 = openmc.Cell(fill=u1)
     c4 = openmc.Cell(fill=u2)
-    u3 = openmc.Universe(cells=[c3, c4])
+    u3 = openmc.DAGMCUniverse(filename="")
+    c5 = openmc.Cell(fill=u3)
+    u4 = openmc.Universe(cells=[c3, c4, c5])
 
-    univs = set(u3.get_all_universes().values())
-    assert not (univs ^ {u1, u2})
+    univs = set(u4.get_all_universes().values())
+    assert not (univs ^ {u1, u2, u3})
 
 
 def test_clone():
@@ -107,11 +109,13 @@ def test_clone():
     c2.fill = openmc.Material()
     c3 = openmc.Cell()
     u1 = openmc.Universe(name='cool', cells=(c1, c2, c3))
+    u1.volume = 1.
 
     u2 = u1.clone()
     assert u2.name == u1.name
     assert u2.cells != u1.cells
     assert u2.get_all_materials() != u1.get_all_materials()
+    assert u2.volume == u1.volume
 
     u2 = u1.clone(clone_materials=False)
     assert u2.get_all_materials() == u1.get_all_materials()
@@ -120,6 +124,33 @@ def test_clone():
     assert next(iter(u3.cells.values())).region ==\
          next(iter(u1.cells.values())).region
 
+    # Change attributes, make sure clone stays intact
+    u1.volume = 2.
+    u1.name = "different name"
+    assert u3.volume != u1.volume
+    assert u3.name != u1.name
+
+    # Test cloning a DAGMC universe
+    dagmc_u = openmc.DAGMCUniverse(filename="", name="DAGMC universe")
+    dagmc_u.volume = 1.
+    dagmc_u.auto_geom_ids = True
+    dagmc_u.auto_mat_ids = True
+    dagmc_u1 = dagmc_u.clone()
+    assert dagmc_u1.name == dagmc_u.name
+    assert dagmc_u1.volume == dagmc_u.volume
+    assert dagmc_u1.auto_geom_ids == dagmc_u.auto_geom_ids
+    assert dagmc_u1.auto_mat_ids == dagmc_u.auto_mat_ids
+
+    # Change attributes, check the clone remained intact
+    dagmc_u.name = "another name"
+    dagmc_u.auto_geom_ids = False
+    dagmc_u.auto_mat_ids = False
+    dagmc_u.volume = 2.
+    assert dagmc_u1.name != dagmc_u.name
+    assert dagmc_u1.volume != dagmc_u.volume
+    assert dagmc_u1.auto_geom_ids != dagmc_u.auto_geom_ids
+    assert dagmc_u1.auto_mat_ids != dagmc_u.auto_mat_ids
+
 
 def test_create_xml(cell_with_lattice):
     cells = [openmc.Cell() for i in range(5)]

tests/unit_tests/weightwindows/test.py

@@ -225,7 +225,7 @@ def test_lower_ww_bounds_shape():
     assert ww.lower_ww_bounds.shape == (2, 3, 4, 1)
 
 
-def test_roundtrip(model, wws):
+def test_roundtrip(run_in_tmpdir, model, wws):
     model.settings.weight_windows = wws
 
     # write the model with weight windows to XML