Add many tests and fix a few bugs (#48)

Author: paulromano

src/openmc_mcnp_adapter/openmc_conversion.py

@@ -334,7 +334,8 @@ def flip_sense(surf):
                     # decide if we want the up or down part of the
                     # cone since one sheet is used
                     up = grad >= 0
-                    surf = cls_cone(z0=offset, r2=angle, up=up)
+                    kwargs = {f"{s['mnemonic']}0": offset, "r2": angle, "up": up}
+                    surf = cls_cone(**kwargs)
             else:
                 raise NotImplementedError(f"{s['mnemonic']} surface with {len(coeffs)} parameters")
         elif s['mnemonic'] == 'rcc':
@@ -919,13 +920,13 @@ def mcnp_to_model(filename, merge_surfaces: bool = True, expand_elements: bool =
     return openmc.Model(geometry, materials, settings)
 
 
-def mcnp_str_to_model(text: str, merge_surfaces: bool = True):
+def mcnp_str_to_model(text: str, **kwargs):
     # Write string to a temporary file
     with tempfile.NamedTemporaryFile('w', delete=False) as fp:
         fp.write(text)
 
     # Parse model from file
-    model = mcnp_to_model(fp.name, merge_surfaces)
+    model = mcnp_to_model(fp.name, **kwargs)
 
     # Remove temporary file and return model
     os.remove(fp.name)

src/openmc_mcnp_adapter/parse.py

@@ -13,7 +13,7 @@
 _KEYWORDS = [
     r'\*?trcl', r'\*?fill', 'tmp', 'u', 'lat',
     'imp:.', 'vol', 'pwt', 'ext:.', 'fcl', 'wwn', 'dxc', 'nonu', 'pd',
-    'elpt', 'cosy', 'bflcl', 'unc',
+    'elpt', 'cosy', 'bflcl', 'unc', 'mat', 'rho',
     'pmt'  # D1SUNED-specific
 ]
 _ANY_KEYWORD = '|'.join(f'(?:{k})' for k in _KEYWORDS)
@@ -44,9 +44,29 @@
 _SAB_RE = re.compile(r'\s*[Mm][Tt](\d+)((?:\s+\S+)+)')
 _MODE_RE = re.compile(r'\s*mode(?:\s+\S+)*')
 _COMPLEMENT_RE = re.compile(r'(#)[ ]*(\d+)')
-_REPEAT_RE = re.compile(r'(\d+)\s+(\d+)[rR]')
 _NUM_RE = re.compile(r'(\d)([+-])(\d)')
 
+_HAS_REPEAT_RE = re.compile(r'\b\d+[rR]\b')
+
+_REPEAT_RE = re.compile(r"""
+    (?P<value>                  # The numeric value to be repeated
+      [+-]?                     #   Optional sign
+      (?:                       #   Mantissa
+        \d+(?:\.\d*)?           #     Digits with optional fractional part (e.g., 3 or 3. or 3.0)
+        |                       #     or
+        \.\d+                   #     Leading-dot form (e.g., .25)
+      )
+      (?:                       #   Optional exponent
+        [eEdD][+-]?\d+          #     E/D exponent with optional sign (e.g., 1e-3, 2D+3)
+        |                       #     or MCNP "bare" exponent without E/D
+        [+-]\d+                 #     appended sign+digits (e.g., 1.0-3 -> 1.0e-3)
+      )?
+    )
+    \s+                         # One or more spaces between value and count
+    (?P<count>\d+)              # The repeat count
+    [rR]                        # The 'R' or 'r' suffix
+""", re.VERBOSE)
+
 
 def float_(val):
     """Convert scientific notation literals that don't have an 'e' in them to float"""
@@ -346,11 +366,12 @@ def sanitize(section: str) -> str:
     section = re.sub('\n {5}', ' ', section)
 
     # Expand repeated numbers
-    m = _REPEAT_RE.search(section)
-    while m is not None:
-        section = _REPEAT_RE.sub(' '.join((int(m.group(2)) + 1)*[m.group(1)]),
-                                 section, 1)
-        m = _REPEAT_RE.search(section)
+    if _HAS_REPEAT_RE.search(section):
+        section = _REPEAT_RE.sub(
+            lambda m: ' '.join([m.group('value')] * (int(m.group('count')) + 1)),
+            section,
+        )
+
     return section
 
 

tests/test_geometry.py

@@ -0,0 +1,157 @@
+from textwrap import dedent
+
+from pytest import mark, approx, param
+from openmc_mcnp_adapter import mcnp_str_to_model
+
+
+@mark.parametrize("whitespace", ["", " ", "\t"])
+def test_cell_complement(whitespace):
+    # Cell 2 corresponds to r < 2 intersected with z > 0
+    mcnp_str = dedent(f"""
+    title
+    100  1 1.0  +1 : -2
+    2    1 1.0  #{whitespace}100
+
+    1  so 2.0
+    2  pz 0.0
+
+    m1   1001.80c  1.0
+    """)
+    model = mcnp_str_to_model(mcnp_str)
+    cell = model.geometry.get_all_cells()[2]
+
+    # Check various points
+    assert (0., 0., 0.1) in cell.region
+    assert (0., 0., -0.1) not in cell.region
+    assert (0., 0., 1.99) in cell.region
+    assert (0., 0., 2.01) not in cell.region
+    assert (1., 1., 1.) in cell.region
+    assert (2., 0., 1.) not in cell.region
+
+
+def test_likenbut():
+    mcnp_str = dedent("""
+    title
+    1   1 -1.0  -1
+    2   LIKE 1 BUT MAT=2 RHO=-2.0 TRCL=(2.0 0.0 0.0)
+
+    1   so 1.0
+
+    m1   1001.80c  1.0
+    m2   1002.80c  1.0
+    """)
+    model = mcnp_str_to_model(mcnp_str)
+    cell = model.geometry.get_all_cells()[2]
+
+    # Material should be changed to m2
+    mat = cell.fill
+    assert 'H2' in mat.get_nuclide_densities()
+
+    # Density should be 2.0 g/cm3
+    assert mat.get_mass_density() == approx(2.0)
+
+    # Points should correspond to sphere of r=1 centered at (2, 0, 0)
+    assert (2.0, 0.0, 0.0) in cell.region
+    assert (0.0, 0.0, 0.0) not in cell.region
+    assert (2.0, 0.9, 0.0) in cell.region
+    assert (2.0, 1.1, 0.0) not in cell.region
+
+
+@mark.parametrize(
+    "cell_card, surface_cards, points_inside, points_outside",
+    [
+        (
+            "1   1 -1.0  -1 TRCL=(2.0 0.0 0.0)",
+            ("1   so 1.0",),
+            [
+                (2.0, 0.0, 0.0),
+                (2.0, 0.9, 0.0),
+                (2.0, 0.0, 0.9),
+            ],
+            [
+                (0.9, 0.0, 0.0),
+                (2.0, 1.1, 0.0),
+                (2.0, 0.0, 1.1),
+            ],
+        ),
+        (
+            "1 0 -1 TRCL=(1.0 0.0 0.0 0.0 -1.0 0.0 1.0 0.0 0.0 0.0 0.0 1.0)",
+            ("1 rpp -0.5 0.5 -0.25 0.25 -0.1 0.1",),
+            [
+                (1.0, 0.0, 0.0),
+                (1.2, 0.0, 0.0),
+                (1.0, 0.4, 0.0),
+                (1.0, -0.4, 0.0),
+            ],
+            [
+                (0.7, 0.0, 0.0),
+                (1.3, 0.0, 0.0),
+                (1.0, 0.6, 0.0),
+                (1.0, 0.0, 0.2),
+            ],
+        ),
+        (
+            "1 0 -1 *TRCL=(1.0 0.0 0.0 90.0 180.0 90.0 0.0 90.0 90.0 90.0 90.0 0.0)",
+            ("1 rpp -0.5 0.5 -0.25 0.25 -0.1 0.1",),
+            [
+                (1.0, 0.0, 0.0),
+                (1.2, 0.0, 0.0),
+                (1.0, 0.4, 0.0),
+                (1.0, -0.4, 0.0),
+            ],
+            [
+                (0.7, 0.0, 0.0),
+                (1.3, 0.0, 0.0),
+                (1.0, 0.6, 0.0),
+                (1.0, 0.0, 0.2),
+            ],
+        ),
+        param(
+            "1 0 -1 TRCL=(0.0 0.0 0.0 1.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 1.0 -1)",
+            ("1 rpp -0.5 0.5 -0.25 0.25 -0.1 0.1",),
+            [],
+            [],
+            marks=mark.xfail(reason="13-parameter TRCL not yet supported"),
+        ),
+    ],
+)
+def test_trcl(cell_card, surface_cards, points_inside, points_outside):
+    surface_block = "\n".join(surface_cards)
+    mcnp_str = dedent(f"""
+    title
+    {cell_card}
+
+    {surface_block}
+
+    m1   1001.80c  1.0
+    """)
+    model = mcnp_str_to_model(mcnp_str)
+    cell = model.geometry.get_all_cells()[1]
+
+    for point in points_inside:
+        assert point in cell.region
+
+    for point in points_outside:
+        assert point not in cell.region
+
+
+def test_trcl_fill():
+    mcnp_str = dedent("""
+    title
+    1 0 -1 FILL=10 TRCL=(2.0 0.0 0.0)
+    2 0 -2 U=10
+    3 0 +2 U=10
+
+    1 so 10.0
+    2 so 1.0
+
+    m1   1001.80c  1.0
+    """)
+    model = mcnp_str_to_model(mcnp_str)
+    geometry = model.geometry
+    cells = geometry.get_all_cells()
+
+    # Make sure that the cells in universe 10 were shifted
+    assert geometry.find((2.0, 0.0, 0.0))[-1] is cells[2]
+    assert geometry.find((4.0, 0.0, 0.0))[-1] is cells[3]
+    assert geometry.find((0.0, 0.0, 0.0))[-1] is cells[3]

tests/test_material.py

@@ -1,9 +1,24 @@
 import textwrap
 
+import openmc
 from openmc_mcnp_adapter import mcnp_str_to_model
 from pytest import approx
 
 
+def convert_material(mat_card: str, density: float, thermal_card: str = "", **kwargs) -> openmc.Material:
+    mcnp_model = textwrap.dedent(f"""
+        title
+        1   1 {density}    1
+
+        1   px 0.0
+
+        {mat_card}
+        {thermal_card}
+    """)
+    model = mcnp_str_to_model(mcnp_model, **kwargs)
+    return model.materials[0]
+
+
 def test_material_clones():
     mcnp_model = textwrap.dedent("""
         title
@@ -24,3 +39,73 @@ def test_material_clones():
     assert cells[2].fill.id == cells[3].fill.id != 1
     assert cells[1].fill.get_mass_density() == approx(1.0)
     assert cells[2].fill.get_mass_density() == approx(2.0)
+
+
+def test_material_suffixes():
+    # H1 with XS suffix; O16 without suffix
+    mat_card = "m1   1001.80c 2.0  8016 1.0"
+    thermal_card = "mt1  lwtr grph.10t"
+    m = convert_material(mat_card, -2.0, thermal_card)
+    nd = m.get_nuclide_densities()
+    assert set(nd.keys()) == {'H1', 'O16'}
+    assert nd['H1'].percent == approx(2.0)
+    assert nd['H1'].percent_type == 'ao'
+    assert nd['O16'].percent == approx(1.0)
+    assert nd['O16'].percent_type == 'ao'
+
+    # Check S(a,b) tables mapped via get_thermal_name
+    # Access private field because there is no public accessor
+    sab_names = {name for (name, _) in getattr(m, '_sab', [])}
+    assert 'c_H_in_H2O' in sab_names
+    assert 'c_Graphite' in sab_names
+
+
+def test_weight_fractions():
+    # 6000 -> natural C, 5010/5011 -> B-10/B-11; negative => weight
+    mat_card = "m1   6000 -0.12  5010 -0.2  5011 -0.8"
+    m = convert_material(mat_card, -1.0)
+    nd = m.get_nuclide_densities()
+
+    # Natural carbon represented as C0 in OpenMC
+    assert 'C0' in nd and nd['C0'].percent_type == 'wo' and nd['C0'].percent == approx(0.12)
+    assert 'B10' in nd and nd['B10'].percent_type == 'wo' and nd['B10'].percent == approx(0.2)
+    assert 'B11' in nd and nd['B11'].percent_type == 'wo' and nd['B11'].percent == approx(0.8)
+
+
+def test_no_expand_elements():
+    # With expand_elements=False, natural oxygen should be added as O0 directly
+    mat_card = "m1   47000  1.0"
+    m = convert_material(mat_card, -1.0, expand_elements=False)
+    nd = m.get_nuclide_densities()
+    assert 'Ag0' in nd
+    assert nd['Ag0'].percent_type == 'ao'
+    assert nd['Ag0'].percent == approx(1.0)
+
+
+def test_mass_density():
+    mat_card = "m1   3006.80c 0.5 3007.80c 0.5"
+    m = convert_material(mat_card, -3.0)
+    assert m.get_mass_density() == approx(3.0)
+
+
+def test_atom_density():
+    mat_card = "m1   3006.80c 0.5 3007.80c 0.5"
+    m = convert_material(mat_card, 0.02)
+    nuclide_densities = m.get_nuclide_atom_densities()
+    assert sum(nuclide_densities.values()) == approx(0.02)
+    assert nuclide_densities['Li6'] == approx(0.01)
+    assert nuclide_densities['Li7'] == approx(0.01)
+
+
+def test_density_no_whitespace():
+    mcnp_model = textwrap.dedent("""
+        title
+        1  1 -4.5( -1 )
+
+        1  px 0.0
+
+        m1  1001.80c  2.0
+    """)
+    model = mcnp_str_to_model(mcnp_model)
+    m = model.materials[0]
+    assert m.get_mass_density() == approx(4.5)