diff --git a/src/openmc_mcnp_adapter/openmc_conversion.py b/src/openmc_mcnp_adapter/openmc_conversion.py index 40deec5..69b2bdf 100644 --- a/src/openmc_mcnp_adapter/openmc_conversion.py +++ b/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) diff --git a/src/openmc_mcnp_adapter/parse.py b/src/openmc_mcnp_adapter/parse.py index 8bf72ad..05b8561 100644 --- a/src/openmc_mcnp_adapter/parse.py +++ b/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 # 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\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 diff --git a/tests/test_geometry.py b/tests/test_geometry.py new file mode 100644 index 0000000..9de30ed --- /dev/null +++ b/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] diff --git a/tests/test_material.py b/tests/test_material.py index 4d7d4da..12dff0a 100644 --- a/tests/test_material.py +++ b/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) diff --git a/tests/test_surfaces.py b/tests/test_surfaces.py index c2909c9..b96fa92 100644 --- a/tests/test_surfaces.py +++ b/tests/test_surfaces.py @@ -1,9 +1,11 @@ from collections.abc import Sequence from textwrap import dedent +import numpy as np import openmc from openmc.model.surface_composite import OrthogonalBox, \ - RectangularParallelepiped, RightCircularCylinder, ConicalFrustum + RectangularParallelepiped, RightCircularCylinder, ConicalFrustum, \ + XConeOneSided, YConeOneSided, ZConeOneSided from openmc_mcnp_adapter import mcnp_str_to_model, get_openmc_surfaces from pytest import approx, mark @@ -28,6 +30,20 @@ def convert_surface(mnemonic: str, params: Sequence[float]) -> openmc.Surface: return surfaces[1] +def test_reflective_surface(): + mcnp_str = dedent(""" + title + 1 1.0 -1 + + *1 so 2.0 + + m1 1001.80c 1.0 + """) + model = mcnp_str_to_model(mcnp_str) + surf = model.geometry.get_all_surfaces()[1] + assert surf.boundary_type == 'reflective' + + @mark.parametrize( "mnemonic, params, expected_type, attrs", [ @@ -44,7 +60,7 @@ def test_planes(mnemonic, params, expected_type, attrs): assert getattr(surf, attr) == approx(value) -def test_plane_from_points(): +def test_plane_9points(): # Points defining plane y = x - 1 coeffs = (1.0, 0.0, 0.0, 2.0, 1.0, 0.0, @@ -57,12 +73,79 @@ def test_plane_from_points(): assert surf.d == approx(1.0) +def test_plane_sense_rule1(): + # In general, origin is required to have negative sense + coeffs = ( + 0., 1., 0., + 1., 1., 0., + 0., 1., 1. + ) + surf = convert_surface("p", coeffs) + assert (0., 0., 0.) in -surf + + +def test_plane_sense_rule2(): + # If plane passes through origin, the point (0, 0, ∞) has positive sense. In + # this case, we use the surface x - z = 0 + coeffs = ( + 1., 0., 1., + 0., 0., 0., + 1., 1., 1. + ) + surf = convert_surface("p", coeffs) + assert (0., 0., 1e10) in +surf + + +def test_plane_sense_rule3(): + # If D = C = 0, the point (0, ∞, 0) has positive sense. In this case, we use + # the surface, x - y = 0 + coeffs = ( + 0., 0., 0., + 1., 1., 0., + 0., 0., 1. + ) + surf = convert_surface("p", coeffs) + assert (0., 1e10, 0.) in +surf + + +def test_plane_sense_rule4(): + # If D = C = B = 0, the point (∞, 0, 0) has positive sense. In this case, we + # use the surface x = 0 + coeffs = ( + 0., 1., 1., + 0., 0., 0., + 0., 0., 1. + ) + surf = convert_surface("p", coeffs) + assert (1e10, 0., 0.) in +surf + + +def test_surface_transformation_with_tr_card(): + mcnp_str = dedent(""" + title + 1 0 -1 + + 1 1 pz 0.0 + + tr1 0.0 1.0 1.0 1.0 0.0 0.0 0.0 0.0 1.0 0.0 -1.0 0.0 + """) + model = mcnp_str_to_model(mcnp_str) + surf = model.geometry.get_all_surfaces()[1] + + # The transformed plane corresponds to -y + 1 = 0, so y > 1 is negative sense + assert (0.0, 1.1, 1.0) in -surf + assert (0.0, 0.9, 1.0) in +surf + + @mark.parametrize( "mnemonic, params", [ ("so", (2.4,)), ("s", (1.0, -2.0, 3.0, 4.5)), ("sph", (1.0, -2.0, 3.0, 4.5)), + ("sx", (2.0, 1.5)), + ("sy", (9.0, 0.5)), + ("sz", (-4.0, 4.0)), ], ) def test_spheres(mnemonic, params): @@ -75,13 +158,28 @@ def test_spheres(mnemonic, params): assert surf.x0 == 0.0 assert surf.y0 == 0.0 assert surf.z0 == 0.0 - assert surf.r == approx(r) + elif mnemonic in ("sx", "sy", "sz"): + center, r = params + if mnemonic == "sx": + assert surf.x0 == approx(center) + assert surf.y0 == 0.0 + assert surf.z0 == 0.0 + elif mnemonic == "sy": + assert surf.x0 == 0.0 + assert surf.y0 == approx(center) + assert surf.z0 == 0.0 + else: + assert surf.x0 == 0.0 + assert surf.y0 == 0.0 + assert surf.z0 == approx(center) else: x0, y0, z0, r = params assert surf.x0 == approx(x0) assert surf.y0 == approx(y0) assert surf.z0 == approx(z0) - assert surf.r == approx(r) + + # Check radius + assert surf.r == approx(r) @mark.parametrize( @@ -209,6 +307,83 @@ def test_axisymmetric_surfaces(mnemonic, params, expected_type, attr, value): assert getattr(surf, attr) == approx(value) +@mark.parametrize("mnemonic", ["x", "y", "z"]) +def test_axisymmetric_surfaces_cone(mnemonic): + # cone with a r=1 bottom at plane=0 and a r=2 top at plane=3 + coeffs = (0.0, 1.0, 3.0, 2.0) + surf = convert_surface(mnemonic, coeffs) + + # Helper to build a point (x,y,z) given radial distance r and axial coord a + def pt(r: float, a: float): + if mnemonic == "x": + return (a, r, 0.0) # axial along x; radius in y + elif mnemonic == "y": + return (r, a, 0.0) # axial along y; radius in x + else: # "z" + return (r, 0.0, a) # axial along z; radius in x + + # Points near the r=1 slice (at axial ~ 0) + assert pt(0.0, 0.01) in -surf + assert pt(0.0, -0.01) in -surf + assert pt(0.99, 0.01) in -surf + assert pt(1.05, 0.01) in +surf + + # Points near the r=2 slice (at axial ~ 3) + assert pt(1.99, 2.99) in -surf + assert pt(2.0, 2.99) in +surf + assert pt(0.0, 2.99) in -surf + assert pt(0.0, 3.01) in -surf + + # Points between the r=1 and r=2 slices (at axial = 1.5) + assert pt(1.49, 1.5) in -surf + assert pt(1.51, 1.5) in +surf + + +@mark.parametrize( + "mnemonic, params, expected_type, up_expected", + [ + # k/x with one-sided flag: +1 -> up=True, -1 -> up=False + ("k/x", (0.0, 0.0, 0.0, 1.0, 1.0), XConeOneSided, True), + ("k/x", (0.0, 0.0, 0.0, 2.5, -1.0), XConeOneSided, False), + ("k/y", (1.0, -2.0, 3.0, 0.5, 1.0), YConeOneSided, True), + ("k/z", (-1.0, 2.0, -3.0, 4.0, -2.0), ZConeOneSided, False), + ], +) +def test_cones_one_sided(mnemonic, params, expected_type, up_expected): + surf = convert_surface(mnemonic, params) + assert isinstance(surf, expected_type) + assert getattr(surf, "up") is up_expected + # Access nested cone attributes for center and r2 + assert surf.cone.x0 == approx(params[0]) + assert surf.cone.y0 == approx(params[1]) + assert surf.cone.z0 == approx(params[2]) + assert surf.cone.r2 == approx(params[3]) + + +@mark.parametrize( + "mnemonic, params, expected_type, up_expected", + [ + ("kx", (0.5, 2.0, 1.0), XConeOneSided, True), + ("kx", (0.5, 2.0, -1.0), XConeOneSided, False), + ("ky", (-0.3, 1.2, 2.0), YConeOneSided, True), + ("ky", (-0.3, 1.2, -2.0), YConeOneSided, False), + ("kz", (2.2, 0.7, 3.0), ZConeOneSided, True), + ("kz", (2.2, 0.7, -3.0), ZConeOneSided, False), + ], +) +def test_cones_one_sided_short(mnemonic, params, expected_type, up_expected): + surf = convert_surface(mnemonic, params) + assert isinstance(surf, expected_type) + assert getattr(surf, "up") is up_expected + + # Check plane position coincides with axis center value + axis = mnemonic[1] # 'x', 'y', or 'z' + assert getattr(surf.plane, f"{axis}0") == approx(params[0]) + + # Check cone radius-squared parameter + assert surf.cone.r2 == approx(params[1]) + + def test_box_macrobody(): coeffs = (0.0, 0.0, 0.0, 1.0, 0.0, 0.0, @@ -259,6 +434,27 @@ def test_rcc_macrobody(coeffs, expected_bottom, expected_top, r, coeff): assert surf.top.d / getattr(surf.top, coeff) == approx(expected_top) +def test_rcc_macrobody_non_axis_aligned(): + coeffs = (0.0, 0.0, 0.0, 1.0, 1.0, 0.0, 0.5) + surf = convert_surface("rcc", coeffs) + assert isinstance(surf, RightCircularCylinder) + + # Points along and around the axis defined by the height vector (1, 1, 0) + midpoint = (0.5, 0.5, 0.0) + axis_unit = np.array(coeffs[3:6]) + + inside_point = (midpoint[0], midpoint[1], 0.1) + outside_radial = (midpoint[0], midpoint[1], 0.6) + below_bottom = -0.1 * axis_unit + top = np.array((1.0, 1.0, 0.0)) + above_top = top + 0.1 * axis_unit + + assert inside_point in -surf + assert outside_radial in +surf + assert below_bottom in +surf + assert above_top in +surf + + def test_trc_macrobody(): coeffs = (0.0, 0.0, 0.0, 0.0, 0.0, 10.0, 1.0, 2.0) surf = convert_surface("trc", coeffs) diff --git a/tests/test_syntax.py b/tests/test_syntax.py new file mode 100644 index 0000000..9dca129 --- /dev/null +++ b/tests/test_syntax.py @@ -0,0 +1,51 @@ +from textwrap import dedent + +from pytest import approx +from openmc_mcnp_adapter import mcnp_str_to_model + + +def test_repeat_shortcut(): + mcnp_str = dedent(""" + title + 1 0 -1 + + 1 gq 1.0 3r 0.0 5r + + m1 1001.80c 3.0 + """) + model = mcnp_str_to_model(mcnp_str) + surf = model.geometry.get_all_surfaces()[1] + print(surf) + + # Make sure A, B, C, and D parameters are 1.0 + for attr in 'abcd': + assert getattr(surf, attr) == 1.0 + + # Make sure E, F, G, H, J, and K parameters are 0.0 + for attr in 'efghj': + assert getattr(surf, attr) == 0.0 + + +def test_comments(): + mcnp_str = dedent(""" + title + c This is a comment line + 1 1 -5.0 -1 $ This is an end-of-line comment + + c surface block + 1 so 3.0 $ Another comment + + m1 1001.80c 3.0 $ Material comment + """) + model = mcnp_str_to_model(mcnp_str) + surf = model.geometry.get_all_surfaces()[1] + cell = model.geometry.get_all_cells()[1] + mat = model.materials[0] + + # Sanity checks + assert surf.r == 3.0 + assert surf.x0 == 0.0 + assert cell.id == 1 + assert str(cell.region) == "-1" + assert 'H1' in mat.get_nuclide_densities() + assert mat.get_mass_density() == approx(5.0)