From 105ea725b25f55e61a59afbaef6098fd34b1f031 Mon Sep 17 00:00:00 2001 From: Massimo Larsen Date: Thu, 14 May 2026 11:15:04 -0700 Subject: [PATCH 01/15] Added bounded surface tallies. Added unit tests for multiple different tallies --- mcdc/numba_types.py | 10 ++ mcdc/object_/tally.py | 80 ++++++++++++++ mcdc/transport/geometry/interface.py | 42 ++++++++ test/unit/transport/tally/__init__.py | 0 test/unit/transport/tally/conftest.py | 100 ++++++++++++++++++ .../transport/tally/test_collision_tally.py | 26 +++++ .../tally/test_nonbounded_surface_tally.py | 54 ++++++++++ .../tally/test_surface_bounds_edge.py | 48 +++++++++ .../tally/test_surface_bounds_validation.py | 33 ++++++ .../tally/test_tally_factory_routing.py | 30 ++++++ .../transport/tally/test_tracklength_tally.py | 28 +++++ 11 files changed, 451 insertions(+) create mode 100644 test/unit/transport/tally/__init__.py create mode 100644 test/unit/transport/tally/conftest.py create mode 100644 test/unit/transport/tally/test_collision_tally.py create mode 100644 test/unit/transport/tally/test_nonbounded_surface_tally.py create mode 100644 test/unit/transport/tally/test_surface_bounds_edge.py create mode 100644 test/unit/transport/tally/test_surface_bounds_validation.py create mode 100644 test/unit/transport/tally/test_tally_factory_routing.py create mode 100644 test/unit/transport/tally/test_tracklength_tally.py diff --git a/mcdc/numba_types.py b/mcdc/numba_types.py index 004a10908..1cc939d25 100644 --- a/mcdc/numba_types.py +++ b/mcdc/numba_types.py @@ -675,6 +675,16 @@ surface_tally = into_dtype([ ('surface_ID', int64), + ('filter_surface_bounds', bool), + ('has_x_bounds', bool), + ('has_y_bounds', bool), + ('has_z_bounds', bool), + ('x_min', float64), + ('x_max', float64), + ('y_min', float64), + ('y_max', float64), + ('z_min', float64), + ('z_max', float64), ('ID', int64), ('parent_ID', int64), ]) diff --git a/mcdc/object_/tally.py b/mcdc/object_/tally.py index f55632bbd..993a0b2db 100644 --- a/mcdc/object_/tally.py +++ b/mcdc/object_/tally.py @@ -35,6 +35,9 @@ SPATIAL_FILTER_CELL, SPATIAL_FILTER_MESH, SPATIAL_FILTER_NONE, + SURFACE_PLANE_X, + SURFACE_PLANE_Y, + SURFACE_PLANE_Z, TALLY_SURFACE, TALLY_COLLISION, TALLY_TRACKLENGTH, @@ -91,6 +94,9 @@ def __new__( polar_reference: Iterable[float] | NoneType = None, energy: Iterable[float] | str | NoneType = None, time: Iterable[float] | NoneType = None, + x: Iterable[float] | NoneType = None, + y: Iterable[float] | NoneType = None, + z: Iterable[float] | NoneType = None, ) -> TallySurface | TallyTracklength | TallyCollision: # Determine type and create the tally self based on the provided # spatial filters and scores @@ -137,6 +143,9 @@ def __init__( polar_reference: Iterable[float] | NoneType = None, energy: Iterable[float] | str | NoneType = None, time: Iterable[float] | NoneType = None, + x: Iterable[float] | NoneType = None, + y: Iterable[float] | NoneType = None, + z: Iterable[float] | NoneType = None, spatial_shape: tuple[int] | NoneType = None, ): # Set name @@ -299,6 +308,16 @@ class TallySurface(Tally): label: str = "surface_tally" # surface: Surface + filter_surface_bounds: bool + has_x_bounds: bool + has_y_bounds: bool + has_z_bounds: bool + x_min: float + x_max: float + y_min: float + y_max: float + z_min: float + z_max: float def __init__( self, @@ -310,6 +329,9 @@ def __init__( polar_reference: Iterable[float] | NoneType = None, energy: Iterable[float] | str | NoneType = None, time: Iterable[float] | NoneType = None, + x: Iterable[float] | NoneType = None, + y: Iterable[float] | NoneType = None, + z: Iterable[float] | NoneType = None, ): type_ = TALLY_SURFACE super(Tally, self).__init__(type_) @@ -333,6 +355,64 @@ def __init__( self.surface = surface surface.tallies.append(self) + # Optional bounds for axis-aligned planar surface tallies. + self.filter_surface_bounds = False + self.has_x_bounds = False + self.has_y_bounds = False + self.has_z_bounds = False + self.x_min = -INF + self.x_max = INF + self.y_min = -INF + self.y_max = INF + self.z_min = -INF + self.z_max = INF + + if x is not None: + x = np.array(x, dtype=float) + if len(x) != 2: + print_error( + "Surface tally bound x must have exactly two values [min, max]." + ) + self.has_x_bounds = True + self.x_min = min(x[0], x[1]) + self.x_max = max(x[0], x[1]) + + if y is not None: + y = np.array(y, dtype=float) + if len(y) != 2: + print_error( + "Surface tally bound y must have exactly two values [min, max]." + ) + self.has_y_bounds = True + self.y_min = min(y[0], y[1]) + self.y_max = max(y[0], y[1]) + + if z is not None: + z = np.array(z, dtype=float) + if len(z) != 2: + print_error( + "Surface tally bound z must have exactly two values [min, max]." + ) + self.has_z_bounds = True + self.z_min = min(z[0], z[1]) + self.z_max = max(z[0], z[1]) + + self.filter_surface_bounds = ( + self.has_x_bounds or self.has_y_bounds or self.has_z_bounds + ) + + if self.filter_surface_bounds: + if surface.type not in (SURFACE_PLANE_X, SURFACE_PLANE_Y, SURFACE_PLANE_Z): + print_error( + "Bounded surface tally currently supports only PlaneX, PlaneY, and PlaneZ surfaces." + ) + if surface.type == SURFACE_PLANE_X and self.has_x_bounds: + print_error("PlaneX surface tally bounds may only use y and/or z.") + if surface.type == SURFACE_PLANE_Y and self.has_y_bounds: + print_error("PlaneY surface tally bounds may only use x and/or z.") + if surface.type == SURFACE_PLANE_Z and self.has_z_bounds: + print_error("PlaneZ surface tally bounds may only use x and/or y.") + def __repr__(self): text = super().__repr__() text += f" - Surface: {self.surface.name}\n" diff --git a/mcdc/transport/geometry/interface.py b/mcdc/transport/geometry/interface.py index 0de3f79aa..05231e060 100644 --- a/mcdc/transport/geometry/interface.py +++ b/mcdc/transport/geometry/interface.py @@ -454,6 +454,13 @@ def surface_crossing(P_arr, simulation, data): for i in range(surface["N_tally"]): tally_ID = int(mcdc_get.surface.tally_IDs(i, surface, data)) tally = simulation["surface_tallies"][tally_ID] + + # Optional bounded surface tally: score only if the crossing point is + # within the specified in-plane bounds. + if tally["filter_surface_bounds"]: + if not _surface_tally_crossing_in_bounds(P, surface, tally): + continue + tally_module.score.surface_tally(P_arr, surface, tally, simulation, data) # Need to check new cell later? @@ -462,6 +469,41 @@ def surface_crossing(P_arr, simulation, data): P["material_ID"] = -1 +@njit +def _surface_tally_crossing_in_bounds(P, surface, tally): + surface_type = surface["type"] + + if surface_type == SURFACE_PLANE_X: + if tally["has_y_bounds"]: + if P["y"] < tally["y_min"] or P["y"] > tally["y_max"]: + return False + if tally["has_z_bounds"]: + if P["z"] < tally["z_min"] or P["z"] > tally["z_max"]: + return False + return True + + if surface_type == SURFACE_PLANE_Y: + if tally["has_x_bounds"]: + if P["x"] < tally["x_min"] or P["x"] > tally["x_max"]: + return False + if tally["has_z_bounds"]: + if P["z"] < tally["z_min"] or P["z"] > tally["z_max"]: + return False + return True + + if surface_type == SURFACE_PLANE_Z: + if tally["has_x_bounds"]: + if P["x"] < tally["x_min"] or P["x"] > tally["x_max"]: + return False + if tally["has_y_bounds"]: + if P["y"] < tally["y_min"] or P["y"] > tally["y_max"]: + return False + return True + + # Constructor validation should prevent non-planar bounded surface tallies. + return True + + # ====================================================================================== # Miscellanies # ====================================================================================== diff --git a/test/unit/transport/tally/__init__.py b/test/unit/transport/tally/__init__.py new file mode 100644 index 000000000..e69de29bb diff --git a/test/unit/transport/tally/conftest.py b/test/unit/transport/tally/conftest.py new file mode 100644 index 000000000..317676b85 --- /dev/null +++ b/test/unit/transport/tally/conftest.py @@ -0,0 +1,100 @@ +import os + +import numpy as np +import pytest + +# Match distribution tests: run Numba-decorated routines in pure Python for unit tests. +os.environ.setdefault("NUMBA_DISABLE_JIT", "1") + +import mcdc +import mcdc.numba_types as type_ +from mcdc.constant import PARTICLE_NEUTRON +from mcdc.main import preparation +from mcdc.object_.simulation import simulation + + +@pytest.fixture(autouse=True) +def reset_simulation(): + # Keep simulation state isolated per test. + simulation.__init__() + yield + simulation.__init__() + + +@pytest.fixture +def material_mg(): + # Minimal multigroup material so particle speed is defined. + return mcdc.MaterialMG(capture=np.array([1.0])) + + +@pytest.fixture +def slab_plane_x(material_mg): + # Two-cell 1D slab with shared interior PlaneX. + s_left = mcdc.Surface.PlaneX(x=-1.0, boundary_condition="vacuum") + s_mid = mcdc.Surface.PlaneX(x=0.0) + s_right = mcdc.Surface.PlaneX(x=1.0, boundary_condition="vacuum") + c_left = mcdc.Cell(region=+s_left & -s_mid, fill=material_mg) + c_right = mcdc.Cell(region=+s_mid & -s_right, fill=material_mg) + return { + "s_left": s_left, + "s_mid": s_mid, + "s_right": s_right, + "c_left": c_left, + "c_right": c_right, + } + + +@pytest.fixture +def surface_tally_context(slab_plane_x): + s_mid = slab_plane_x["s_mid"] + + # Same surface, with and without explicit y/z bounds. + unbounded_tally_obj = mcdc.Tally(surface=s_mid, scores=["net-current"]) + bounded_tally_obj = mcdc.Tally( + surface=s_mid, + y=[-0.5, 0.5], + z=[-0.25, 0.25], + scores=["net-current"], + ) + + # Build compiled structures. + mcdc_container, data = preparation() + mcdc_struct = mcdc_container[0] + + # Compiled tally handles. + unbounded_tally = mcdc_struct["surface_tallies"][unbounded_tally_obj.child_ID] + bounded_tally = mcdc_struct["surface_tallies"][bounded_tally_obj.child_ID] + + # Particle for direct crossing-based tally-kernel testing. + particle_container = np.zeros(1, type_.particle) + particle = particle_container[0] + particle["particle_type"] = PARTICLE_NEUTRON + particle["material_ID"] = 0 + particle["g"] = 0 + particle["x"] = 0.0 + particle["y"] = 0.0 + particle["z"] = 0.0 + particle["t"] = 0.0 + particle["ux"] = 0.5 + particle["uy"] = 0.0 + particle["uz"] = 0.0 + particle["w"] = 2.0 + + return { + "data": data, + "mcdc_struct": mcdc_struct, + "particle": particle, + "particle_container": particle_container, + "s_mid": s_mid, + "unbounded_tally": unbounded_tally, + "bounded_tally": bounded_tally, + } + + +@pytest.fixture +def bin_value(): + def _value(surface_tally, mcdc_struct, data): + tally_base = mcdc_struct["tallies"][surface_tally["parent_ID"]] + return data[tally_base["bin_offset"]] + + return _value diff --git a/test/unit/transport/tally/test_collision_tally.py b/test/unit/transport/tally/test_collision_tally.py new file mode 100644 index 000000000..47927dc01 --- /dev/null +++ b/test/unit/transport/tally/test_collision_tally.py @@ -0,0 +1,26 @@ +import pytest + +import mcdc +from mcdc.constant import SPATIAL_FILTER_MESH +from mcdc.object_.tally import TallyCollision + + +def test_collision_tally_with_mesh_filter(): + mesh = mcdc.MeshUniform( + "mesh", + x=(-1.0, 0.5, 2), + y=(-1.0, 1.0, 1), + z=(-1.0, 1.0, 1), + ) + tally = mcdc.Tally(mesh=mesh, scores=["energy_deposition"]) + + assert isinstance(tally, TallyCollision) + assert tally.spatial_filter_type == SPATIAL_FILTER_MESH + assert tally.spatial_filter_ID == mesh.ID + + +def test_collision_tally_requires_mesh(capsys): + with pytest.raises(SystemExit): + mcdc.Tally(scores=["energy_deposition"]) + out = capsys.readouterr().out + assert "currently only supported with a mesh spatial filter" in out diff --git a/test/unit/transport/tally/test_nonbounded_surface_tally.py b/test/unit/transport/tally/test_nonbounded_surface_tally.py new file mode 100644 index 000000000..c6ba9c855 --- /dev/null +++ b/test/unit/transport/tally/test_nonbounded_surface_tally.py @@ -0,0 +1,54 @@ +import numpy as np +import pytest + +from mcdc.transport.geometry import interface as geometry_interface + + +def test_surface_tally_bounds_fields(surface_tally_context): + unbounded_tally = surface_tally_context["unbounded_tally"] + bounded_tally = surface_tally_context["bounded_tally"] + + assert unbounded_tally["filter_surface_bounds"] == 0 + assert bounded_tally["filter_surface_bounds"] == 1 + assert bounded_tally["has_x_bounds"] == 0 + assert bounded_tally["has_y_bounds"] == 1 + assert bounded_tally["has_z_bounds"] == 1 + + +@pytest.mark.parametrize( + "ux, y, z, expected_unbounded, expected_bounded", + [ + (0.5, 0.0, 0.0, 2.0, 2.0), # inside bounds: left -> right + (-0.5, 0.0, 0.0, -2.0, -2.0), # inside bounds: right -> left + (0.5, 0.8, 0.0, 2.0, 0.0), # outside y bounds + (0.5, 0.0, 0.3, 2.0, 0.0), # outside z bounds + ], + ids=["inside_l2r", "inside_r2l", "outside_y_bounds", "outside_z_bounds"], +) +def test_unbounded_surface_tally_scoring( + surface_tally_context, + bin_value, + ux, + y, + z, + expected_unbounded, + expected_bounded, +): + data = surface_tally_context["data"] + mcdc_struct = surface_tally_context["mcdc_struct"] + particle = surface_tally_context["particle"] + particle_container = surface_tally_context["particle_container"] + s_mid = surface_tally_context["s_mid"] + unbounded_tally = surface_tally_context["unbounded_tally"] + bounded_tally = surface_tally_context["bounded_tally"] + + particle["alive"] = True + particle["surface_ID"] = s_mid.ID + particle["x"] = 0.0 + particle["y"] = y + particle["z"] = z + particle["ux"] = ux + geometry_interface.surface_crossing(particle_container, mcdc_struct, data) + + assert np.isclose(bin_value(unbounded_tally, mcdc_struct, data), expected_unbounded) + assert np.isclose(bin_value(bounded_tally, mcdc_struct, data), expected_bounded) diff --git a/test/unit/transport/tally/test_surface_bounds_edge.py b/test/unit/transport/tally/test_surface_bounds_edge.py new file mode 100644 index 000000000..9c6e96a7a --- /dev/null +++ b/test/unit/transport/tally/test_surface_bounds_edge.py @@ -0,0 +1,48 @@ +import numpy as np +import pytest + +from mcdc.transport.geometry import interface as geometry_interface + + +@pytest.mark.parametrize( + "y, z, expected_bounded", + [ + (-0.5, 0.0, 2.0), # exactly at y lower bound: included + (0.5, 0.0, 2.0), # exactly at y upper bound: included + (-0.500001, 0.0, 0.0), # just outside y lower bound: excluded + (0.500001, 0.0, 0.0), # just outside y upper bound: excluded + (0.0, -0.25, 2.0), # exactly at z lower bound: included + (0.0, 0.25, 2.0), # exactly at z upper bound: included + (0.0, -0.250001, 0.0), # just outside z lower bound: excluded + (0.0, 0.250001, 0.0), # just outside z upper bound: excluded + ], + ids=[ + "y_at_lower", + "y_at_upper", + "y_outside_lower", + "y_outside_upper", + "z_at_lower", + "z_at_upper", + "z_outside_lower", + "z_outside_upper", + ], +) +def test_surface_bounds_edges(surface_tally_context, bin_value, y, z, expected_bounded): + data = surface_tally_context["data"] + mcdc_struct = surface_tally_context["mcdc_struct"] + particle = surface_tally_context["particle"] + particle_container = surface_tally_context["particle_container"] + s_mid = surface_tally_context["s_mid"] + unbounded_tally = surface_tally_context["unbounded_tally"] + bounded_tally = surface_tally_context["bounded_tally"] + + particle["alive"] = True + particle["surface_ID"] = s_mid.ID + particle["x"] = 0.0 + particle["y"] = y + particle["z"] = z + particle["ux"] = 0.5 + geometry_interface.surface_crossing(particle_container, mcdc_struct, data) + + assert np.isclose(bin_value(unbounded_tally, mcdc_struct, data), 2.0) + assert np.isclose(bin_value(bounded_tally, mcdc_struct, data), expected_bounded) diff --git a/test/unit/transport/tally/test_surface_bounds_validation.py b/test/unit/transport/tally/test_surface_bounds_validation.py new file mode 100644 index 000000000..21159e117 --- /dev/null +++ b/test/unit/transport/tally/test_surface_bounds_validation.py @@ -0,0 +1,33 @@ +import pytest + +import mcdc + + +def test_surface_bounds_validation_wrong_axis_on_planex(slab_plane_x, capsys): + with pytest.raises(SystemExit): + mcdc.Tally(surface=slab_plane_x["s_mid"], x=[-1.0, 1.0], scores=["net-current"]) + out = capsys.readouterr().out + assert "PlaneX surface tally bounds may only use y and/or z." in out + + +def test_surface_bounds_validation_bound_shape(slab_plane_x, capsys): + with pytest.raises(SystemExit): + mcdc.Tally( + surface=slab_plane_x["s_mid"], + y=[-1.0, 0.0, 1.0], + scores=["net-current"], + ) + out = capsys.readouterr().out + assert "Surface tally bound y must have exactly two values [min, max]." in out + + +def test_surface_bounds_validation_nonplanar_surface(capsys, material_mg): + s_cyl = mcdc.Surface.CylinderZ(center=(0.0, 0.0), radius=1.0) + mcdc.Cell(region=-s_cyl, fill=material_mg) + with pytest.raises(SystemExit): + mcdc.Tally(surface=s_cyl, x=[-0.5, 0.5], scores=["net-current"]) + out = capsys.readouterr().out + assert ( + "Bounded surface tally currently supports only PlaneX, PlaneY, and PlaneZ" + in out + ) diff --git a/test/unit/transport/tally/test_tally_factory_routing.py b/test/unit/transport/tally/test_tally_factory_routing.py new file mode 100644 index 000000000..51de60daa --- /dev/null +++ b/test/unit/transport/tally/test_tally_factory_routing.py @@ -0,0 +1,30 @@ +import pytest + +import mcdc +from mcdc.constant import SPATIAL_FILTER_NONE +from mcdc.object_.tally import TallyCollision, TallySurface, TallyTracklength + + +def test_tally_factory_routing_surface_vs_tracklength_vs_collision(slab_plane_x): + surface_tally = mcdc.Tally(surface=slab_plane_x["s_mid"], scores=["net-current"]) + tracklength_tally = mcdc.Tally(scores=["flux"]) + + mesh = mcdc.MeshUniform( + "mesh", + x=(-1.0, 0.5, 2), + y=(-1.0, 1.0, 1), + z=(-1.0, 1.0, 1), + ) + collision_tally = mcdc.Tally(mesh=mesh, scores=["energy_deposition"]) + + assert isinstance(surface_tally, TallySurface) + assert isinstance(tracklength_tally, TallyTracklength) + assert isinstance(collision_tally, TallyCollision) + assert tracklength_tally.spatial_filter_type == SPATIAL_FILTER_NONE + + +def test_tally_factory_routing_invalid_score_mix(capsys): + with pytest.raises(SystemExit): + mcdc.Tally(scores=["flux", "energy_deposition"]) + out = capsys.readouterr().out + assert "Cannot mix tracklength scores with collision ones." in out diff --git a/test/unit/transport/tally/test_tracklength_tally.py b/test/unit/transport/tally/test_tracklength_tally.py new file mode 100644 index 000000000..9e61d96fa --- /dev/null +++ b/test/unit/transport/tally/test_tracklength_tally.py @@ -0,0 +1,28 @@ +import mcdc +from mcdc.constant import SPATIAL_FILTER_CELL, SPATIAL_FILTER_MESH +from mcdc.object_.tally import TallyTracklength + + +def test_tracklength_tally_with_cell_filter(slab_plane_x): + tally = mcdc.Tally(cell=slab_plane_x["c_left"], scores=["flux", "capture"]) + + assert isinstance(tally, TallyTracklength) + assert tally.spatial_filter_type == SPATIAL_FILTER_CELL + assert tally.spatial_filter_ID == slab_plane_x["c_left"].ID + + +def test_tracklength_tally_with_mesh_filter(): + mesh = mcdc.MeshUniform( + "mesh", + x=(-1.0, 0.5, 2), + y=(-1.0, 1.0, 1), + z=(-1.0, 1.0, 1), + ) + tally = mcdc.Tally(mesh=mesh, scores=["flux"]) + + assert isinstance(tally, TallyTracklength) + assert tally.spatial_filter_type == SPATIAL_FILTER_MESH + assert tally.spatial_filter_ID == mesh.ID + assert tally.mesh_stride_z == 1 + assert tally.mesh_stride_y == mesh.Nz + assert tally.mesh_stride_x == mesh.Nz * mesh.Ny From f824fdab7970884312f18ff48f1387c6c9481f34 Mon Sep 17 00:00:00 2001 From: Massimo Larsen Date: Sat, 16 May 2026 19:54:15 -0700 Subject: [PATCH 02/15] Add cell current tally --- mcdc/constant.py | 7 +- mcdc/mcdc_get/__init__.py | 2 + mcdc/mcdc_set/__init__.py | 2 + mcdc/numba_types.py | 8 ++ mcdc/object_/tally.py | 117 ++++++++++++++++-- mcdc/transport/geometry/interface.py | 59 ++++++++- mcdc/transport/tally/score.py | 60 +++++++++ test/unit/transport/tally/test_cell_tally.py | 102 +++++++++++++++ .../tally/test_tally_factory_routing.py | 32 ++++- 9 files changed, 372 insertions(+), 17 deletions(-) create mode 100644 test/unit/transport/tally/test_cell_tally.py diff --git a/mcdc/constant.py b/mcdc/constant.py index 8adaa4a6e..27d31c8d3 100644 --- a/mcdc/constant.py +++ b/mcdc/constant.py @@ -5,8 +5,9 @@ # Tally types TALLY_SURFACE = 0 -TALLY_COLLISION = 1 -TALLY_TRACKLENGTH = 2 +TALLY_CELL = 1 +TALLY_COLLISION = 2 +TALLY_TRACKLENGTH = 3 # Tally spatial filter types SPATIAL_FILTER_NONE = 0 @@ -33,6 +34,8 @@ SCORE_TIME_MOMENT_MU_SQ = 11 SCORE_SPACE_MOMENT_MU_SQ = 12 SCORE_ENERGY_DEPOSITION = 13 +SCORE_CURRENT_IN = 14 +SCORE_CURRENT_OUT = 15 # Boundary condition BC_NONE = 0 diff --git a/mcdc/mcdc_get/__init__.py b/mcdc/mcdc_get/__init__.py index 31b31afac..03dbf2cde 100644 --- a/mcdc/mcdc_get/__init__.py +++ b/mcdc/mcdc_get/__init__.py @@ -102,6 +102,8 @@ import mcdc.mcdc_get.surface_tally as surface_tally +import mcdc.mcdc_get.cell_tally as cell_tally + import mcdc.mcdc_get.collision_tally as collision_tally import mcdc.mcdc_get.tracklength_tally as tracklength_tally diff --git a/mcdc/mcdc_set/__init__.py b/mcdc/mcdc_set/__init__.py index e91e83a32..cef3b51c6 100644 --- a/mcdc/mcdc_set/__init__.py +++ b/mcdc/mcdc_set/__init__.py @@ -102,6 +102,8 @@ import mcdc.mcdc_set.surface_tally as surface_tally +import mcdc.mcdc_set.cell_tally as cell_tally + import mcdc.mcdc_set.collision_tally as collision_tally import mcdc.mcdc_set.tracklength_tally as tracklength_tally diff --git a/mcdc/numba_types.py b/mcdc/numba_types.py index 1cc939d25..b47342180 100644 --- a/mcdc/numba_types.py +++ b/mcdc/numba_types.py @@ -689,6 +689,12 @@ ('parent_ID', int64), ]) +cell_tally = into_dtype([ + ('cell_ID', int64), + ('ID', int64), + ('parent_ID', int64), +]) + collision_tally = into_dtype([ ('spatial_filter_type', int64), ('spatial_filter_ID', int64), @@ -838,6 +844,8 @@ def set_simulation(N: dict): ('N_tally', int64), ('surface_tallies', surface_tally, (N['surface_tally'])), ('N_surface_tally', int64), + ('cell_tallies', cell_tally, (N['cell_tally'])), + ('N_cell_tally', int64), ('collision_tallies', collision_tally, (N['collision_tally'])), ('N_collision_tally', int64), ('tracklength_tallies', tracklength_tally, (N['tracklength_tally'])), diff --git a/mcdc/object_/tally.py b/mcdc/object_/tally.py index 993a0b2db..087b4e2b9 100644 --- a/mcdc/object_/tally.py +++ b/mcdc/object_/tally.py @@ -32,12 +32,15 @@ SCORE_FISSION, SCORE_NET_CURRENT, SCORE_ENERGY_DEPOSITION, + SCORE_CURRENT_IN, + SCORE_CURRENT_OUT, SPATIAL_FILTER_CELL, SPATIAL_FILTER_MESH, SPATIAL_FILTER_NONE, SURFACE_PLANE_X, SURFACE_PLANE_Y, SURFACE_PLANE_Z, + TALLY_CELL, TALLY_SURFACE, TALLY_COLLISION, TALLY_TRACKLENGTH, @@ -48,6 +51,7 @@ from mcdc.print_ import print_1d_array, print_error SURFACE_SCORES = set(["net-current"]) +CELL_CURRENT_SCORES = set(["net-current", "current-in", "current-out"]) TRACKLENGTH_SCORES = set(["flux", "density", "collision", "capture", "fission"]) COLLISION_SCORES = set(["energy_deposition"]) @@ -97,19 +101,47 @@ def __new__( x: Iterable[float] | NoneType = None, y: Iterable[float] | NoneType = None, z: Iterable[float] | NoneType = None, - ) -> TallySurface | TallyTracklength | TallyCollision: + ) -> TallySurface | TallyCell | TallyTracklength | TallyCollision: # Determine type and create the tally self based on the provided # spatial filters and scores - # Surface tally - if surface is not None: + has_cell_current_score = any(score in CELL_CURRENT_SCORES for score in scores) + + # Surface/cell net-current tally + if has_cell_current_score: for score in scores: - if not score in SURFACE_SCORES: + if score not in CELL_CURRENT_SCORES: + print_error( + "Cannot mix cell-current scores with non-current scores " + f"in one tally. Cell-current scores: {CELL_CURRENT_SCORES}." + ) + + if surface is not None and cell is not None: + print_error( + "Net-current tally must specify exactly one of surface or cell." + ) + + if surface is None and cell is None: + print_error( + "Scoring 'net-current' needs either a surface or a cell tally." + ) + + if surface is not None: + if not set(scores) <= SURFACE_SCORES: print_error( - f"Scoring '{score}' with surface tally is not supported. " - f"Supported surface tally scores: {SURFACE_SCORES}." + "Surface tally currently supports only " + "scores=['net-current']." ) - return super().__new__(TallySurface) + return super().__new__(TallySurface) + return super().__new__(TallyCell) + + # Unsupported score for explicit surface selector + if surface is not None: + for score in scores: + print_error( + f"Scoring '{score}' with surface tally is not supported. " + f"Supported surface tally scores: {SURFACE_SCORES}." + ) # Collision tally if set(scores) <= COLLISION_SCORES: @@ -119,11 +151,6 @@ def __new__( if set(scores) <= TRACKLENGTH_SCORES: return super().__new__(TallyTracklength) - # Error: Missing surface for surface score - for score in scores: - if score in SURFACE_SCORES and surface is None: - print_error(f"Scoring '{score}' needs a surface tally.") - # Error: Unsupported score combination print_error( "Cannot mix tracklength scores with collision ones." @@ -169,6 +196,10 @@ def __init__( self.scores.append(SCORE_FISSION) elif score == "net-current": self.scores.append(SCORE_NET_CURRENT) + elif score == "current-in": + self.scores.append(SCORE_CURRENT_IN) + elif score == "current-out": + self.scores.append(SCORE_CURRENT_OUT) elif score == "energy_deposition": self.scores.append(SCORE_ENERGY_DEPOSITION) else: @@ -277,6 +308,8 @@ def decode_type(type_): return "Tracklength tally" elif type_ == TALLY_SURFACE: return "Surface tally" + elif type_ == TALLY_CELL: + return "Cell tally" elif type_ == TALLY_COLLISION: return "Collision tally" @@ -296,6 +329,66 @@ def decode_score_type(type_, lower_case=False): return "Net current" if not lower_case else "net-current" elif type_ == SCORE_ENERGY_DEPOSITION: return "Energy deposition" if not lower_case else "energy_deposition" + elif type_ == SCORE_CURRENT_IN: + return "Current in" if not lower_case else "current-in" + elif type_ == SCORE_CURRENT_OUT: + return "Current out" if not lower_case else "current-out" + + +# ====================================================================================== +# Cell tally +# ====================================================================================== + + +class TallyCell(Tally): + # Annotations for Numba mode + label: str = "cell_tally" + # + cell: Cell + + def __init__( + self, + cell: Cell, + name: str = "", + scores: list[str] = ["net-current"], + mu: Iterable[float] | NoneType = None, + azi: Iterable[float] | NoneType = None, + polar_reference: Iterable[float] | NoneType = None, + energy: Iterable[float] | str | NoneType = None, + time: Iterable[float] | NoneType = None, + ): + type_ = TALLY_CELL + super(Tally, self).__init__(type_) + super().__init__( + name, + scores, + mu=mu, + azi=azi, + polar_reference=polar_reference, + energy=energy, + time=time, + ) + + if not set(self.scores) <= { + SCORE_NET_CURRENT, + SCORE_CURRENT_IN, + SCORE_CURRENT_OUT, + }: + print_error( + "Cell tally currently supports only " + "scores=['net-current', 'current-in', 'current-out']." + ) + + # Set cell and attach tally to the cell. + self.cell = cell + cell.tallies.append(self) + + def __repr__(self): + text = super().__repr__() + text += f" - Cell: {self.cell.name}\n" + text += super()._phasespace_filter_text() + text += f" - Bin shape [mu, azi, energy, time, score]: {self.bin_shape} \n" + return text # ====================================================================================== diff --git a/mcdc/transport/geometry/interface.py b/mcdc/transport/geometry/interface.py index 05231e060..6268966e0 100644 --- a/mcdc/transport/geometry/interface.py +++ b/mcdc/transport/geometry/interface.py @@ -441,10 +441,17 @@ def distance_to_nearest_surface(particle_container, cell, simulation, data): @njit def surface_crossing(P_arr, simulation, data): P = P_arr[0] + crossed_surface_ID = P["surface_ID"] - # Apply BC - surface = simulation["surfaces"][P["surface_ID"]] + # Determine pre/post top-cell IDs for non-reflective transfers. + pre_cell_ID = -1 + post_cell_ID = -1 + surface = simulation["surfaces"][crossed_surface_ID] BC = surface["boundary_condition"] + if BC != BC_REFLECTIVE: + pre_cell_ID, post_cell_ID = _get_crossing_top_cell_IDs(P_arr, simulation, data) + + # Apply BC if BC == BC_VACUUM: P["alive"] = False elif BC == BC_REFLECTIVE: @@ -463,12 +470,60 @@ def surface_crossing(P_arr, simulation, data): tally_module.score.surface_tally(P_arr, surface, tally, simulation, data) + # Score cell net-current tallies tied to this crossed surface. + if BC != BC_REFLECTIVE: + for i in range(simulation["N_cell_tally"]): + tally = simulation["cell_tallies"][i] + cell = simulation["cells"][tally["cell_ID"]] + if not _cell_has_surface(cell, crossed_surface_ID, data): + continue + tally_module.score.cell_tally( + P_arr, tally, pre_cell_ID, post_cell_ID, simulation, data + ) + # Need to check new cell later? if P["alive"] and not BC == BC_REFLECTIVE: P["cell_ID"] = -1 P["material_ID"] = -1 +@njit +def _get_crossing_top_cell_IDs(particle_container, simulation, data): + P = particle_container[0] + x = P["x"] + y = P["y"] + z = P["z"] + ux = P["ux"] + uy = P["uy"] + uz = P["uz"] + + epsilon = 10.0 * COINCIDENCE_TOLERANCE + + P["x"] = x - epsilon * ux + P["y"] = y - epsilon * uy + P["z"] = z - epsilon * uz + pre_cell_ID = get_cell(particle_container, UNIVERSE_ROOT, simulation, data) + + P["x"] = x + epsilon * ux + P["y"] = y + epsilon * uy + P["z"] = z + epsilon * uz + post_cell_ID = get_cell(particle_container, UNIVERSE_ROOT, simulation, data) + + P["x"] = x + P["y"] = y + P["z"] = z + + return pre_cell_ID, post_cell_ID + + +@njit +def _cell_has_surface(cell, surface_ID, data): + for i in range(cell["N_surface"]): + if int(mcdc_get.cell.surface_IDs(i, cell, data)) == surface_ID: + return True + return False + + @njit def _surface_tally_crossing_in_bounds(P, surface, tally): surface_type = surface["type"] diff --git a/mcdc/transport/tally/score.py b/mcdc/transport/tally/score.py index 924c03386..8a08deca4 100644 --- a/mcdc/transport/tally/score.py +++ b/mcdc/transport/tally/score.py @@ -25,6 +25,8 @@ SCORE_FISSION, SCORE_NET_CURRENT, SCORE_ENERGY_DEPOSITION, + SCORE_CURRENT_IN, + SCORE_CURRENT_OUT, SPATIAL_FILTER_MESH, ) from mcdc.transport.geometry.surface import get_normal_component @@ -75,6 +77,64 @@ def surface_tally(particle_container, surface, tally, simulation, data): util.atomic_add(data, idx_base + i_score, score) +@njit +def cell_tally( + particle_container, + tally, + pre_cell_ID, + post_cell_ID, + simulation, + data, +): + particle = particle_container[0] + tally_base = simulation["tallies"][tally["parent_ID"]] + + # Get filter indices + MG_mode = simulation["settings"]["neutron_multigroup_mode"] + i_mu, i_azi, i_energy, i_time = get_filter_indices( + particle_container, tally_base, data, MG_mode + ) + + # No score if outside non-changing phase-space bins + if i_mu == -1 or i_azi == -1 or i_energy == -1 or i_time == -1: + return + + # Incoming/outgoing event flags + target_cell_ID = tally["cell_ID"] + entered = pre_cell_ID != target_cell_ID and post_cell_ID == target_cell_ID + exited = pre_cell_ID == target_cell_ID and post_cell_ID != target_cell_ID + if not entered and not exited: + return + + # Tally index + idx_base = ( + tally_base["bin_offset"] + + i_mu * tally_base["stride_mu"] + + i_azi * tally_base["stride_azi"] + + i_energy * tally_base["stride_energy"] + + i_time * tally_base["stride_time"] + ) + + # Score + for i_score in range(tally_base["scores_length"]): + score_type = mcdc_get.tally.scores(i_score, tally_base, data) + score = 0.0 + if score_type == SCORE_NET_CURRENT: + if entered: + score = particle["w"] + elif exited: + score = -particle["w"] + elif score_type == SCORE_CURRENT_IN: + if entered: + score = particle["w"] + elif score_type == SCORE_CURRENT_OUT: + if exited: + score = -particle["w"] + + if score != 0.0: + util.atomic_add(data, idx_base + i_score, score) + + # ====================================================================================== # Collision tally # ====================================================================================== diff --git a/test/unit/transport/tally/test_cell_tally.py b/test/unit/transport/tally/test_cell_tally.py new file mode 100644 index 000000000..2726ad383 --- /dev/null +++ b/test/unit/transport/tally/test_cell_tally.py @@ -0,0 +1,102 @@ +import numpy as np + +import mcdc +import mcdc.numba_types as type_ +from mcdc.constant import PARTICLE_NEUTRON +from mcdc.main import preparation +from mcdc.transport.geometry import interface as geometry_interface + + +def _new_particle(surface_ID, x, ux, w=2.0): + particle_container = np.zeros(1, type_.particle) + particle = particle_container[0] + particle["alive"] = True + particle["particle_type"] = PARTICLE_NEUTRON + particle["surface_ID"] = surface_ID + particle["material_ID"] = 0 + particle["g"] = 0 + particle["x"] = x + particle["y"] = 0.0 + particle["z"] = 0.0 + particle["t"] = 0.0 + particle["ux"] = ux + particle["uy"] = 0.0 + particle["uz"] = 0.0 + particle["w"] = w + return particle_container + + +def _bin_value(cell_tally, mcdc_struct, data): + tally_base = mcdc_struct["tallies"][cell_tally["parent_ID"]] + return data[tally_base["bin_offset"]] + + +def _bin_value_score(cell_tally, mcdc_struct, data, score_idx): + tally_base = mcdc_struct["tallies"][cell_tally["parent_ID"]] + return data[tally_base["bin_offset"] + score_idx] + + +def test_cell_tally_sign_for_incoming_and_outgoing(slab_plane_x): + cell_tally_obj = mcdc.Tally(cell=slab_plane_x["c_right"], scores=["net-current"]) + mcdc_container, data = preparation() + mcdc_struct = mcdc_container[0] + cell_tally = mcdc_struct["cell_tallies"][cell_tally_obj.child_ID] + + # Left -> right across the shared interior surface: incoming to c_right (+) + particle_container = _new_particle(slab_plane_x["s_mid"].ID, x=0.0, ux=0.5) + geometry_interface.surface_crossing(particle_container, mcdc_struct, data) + assert np.isclose(_bin_value(cell_tally, mcdc_struct, data), 2.0) + + # Right -> left across the shared interior surface: outgoing from c_right (-) + particle_container = _new_particle(slab_plane_x["s_mid"].ID, x=0.0, ux=-0.5) + geometry_interface.surface_crossing(particle_container, mcdc_struct, data) + assert np.isclose(_bin_value(cell_tally, mcdc_struct, data), 0.0) + + +def test_cell_tally_records_in_and_out_separately(slab_plane_x): + cell_tally_obj = mcdc.Tally( + cell=slab_plane_x["c_right"], + scores=["net-current", "current-in", "current-out"], + ) + mcdc_container, data = preparation() + mcdc_struct = mcdc_container[0] + cell_tally = mcdc_struct["cell_tallies"][cell_tally_obj.child_ID] + + # One incoming and one outgoing crossing. + particle_container = _new_particle(slab_plane_x["s_mid"].ID, x=0.0, ux=0.5) + geometry_interface.surface_crossing(particle_container, mcdc_struct, data) + particle_container = _new_particle(slab_plane_x["s_mid"].ID, x=0.0, ux=-0.5) + geometry_interface.surface_crossing(particle_container, mcdc_struct, data) + + # Scores are in requested order: [net-current, current-in, current-out]. + # current-out uses negative sign for exiting particles. + assert np.isclose(_bin_value_score(cell_tally, mcdc_struct, data, 0), 0.0) + assert np.isclose(_bin_value_score(cell_tally, mcdc_struct, data, 1), 2.0) + assert np.isclose(_bin_value_score(cell_tally, mcdc_struct, data, 2), -2.0) + + +def test_cell_tally_counts_outgoing_to_vacuum(slab_plane_x): + cell_tally_obj = mcdc.Tally(cell=slab_plane_x["c_right"], scores=["net-current"]) + mcdc_container, data = preparation() + mcdc_struct = mcdc_container[0] + cell_tally = mcdc_struct["cell_tallies"][cell_tally_obj.child_ID] + + # c_right -> vacuum across the right boundary: outgoing (-) + particle_container = _new_particle(slab_plane_x["s_right"].ID, x=1.0, ux=0.5) + geometry_interface.surface_crossing(particle_container, mcdc_struct, data) + assert np.isclose(_bin_value(cell_tally, mcdc_struct, data), -2.0) + + +def test_cell_tally_ignores_reflective_crossing(material_mg): + s_left = mcdc.Surface.PlaneX(x=-1.0, boundary_condition="vacuum") + s_right = mcdc.Surface.PlaneX(x=1.0, boundary_condition="reflective") + c_mid = mcdc.Cell(region=+s_left & -s_right, fill=material_mg) + cell_tally_obj = mcdc.Tally(cell=c_mid, scores=["net-current"]) + + mcdc_container, data = preparation() + mcdc_struct = mcdc_container[0] + cell_tally = mcdc_struct["cell_tallies"][cell_tally_obj.child_ID] + + particle_container = _new_particle(s_right.ID, x=1.0, ux=0.5) + geometry_interface.surface_crossing(particle_container, mcdc_struct, data) + assert np.isclose(_bin_value(cell_tally, mcdc_struct, data), 0.0) diff --git a/test/unit/transport/tally/test_tally_factory_routing.py b/test/unit/transport/tally/test_tally_factory_routing.py index 51de60daa..d033b1580 100644 --- a/test/unit/transport/tally/test_tally_factory_routing.py +++ b/test/unit/transport/tally/test_tally_factory_routing.py @@ -2,11 +2,17 @@ import mcdc from mcdc.constant import SPATIAL_FILTER_NONE -from mcdc.object_.tally import TallyCollision, TallySurface, TallyTracklength +from mcdc.object_.tally import ( + TallyCell, + TallyCollision, + TallySurface, + TallyTracklength, +) def test_tally_factory_routing_surface_vs_tracklength_vs_collision(slab_plane_x): surface_tally = mcdc.Tally(surface=slab_plane_x["s_mid"], scores=["net-current"]) + cell_tally = mcdc.Tally(cell=slab_plane_x["c_right"], scores=["net-current"]) tracklength_tally = mcdc.Tally(scores=["flux"]) mesh = mcdc.MeshUniform( @@ -18,6 +24,7 @@ def test_tally_factory_routing_surface_vs_tracklength_vs_collision(slab_plane_x) collision_tally = mcdc.Tally(mesh=mesh, scores=["energy_deposition"]) assert isinstance(surface_tally, TallySurface) + assert isinstance(cell_tally, TallyCell) assert isinstance(tracklength_tally, TallyTracklength) assert isinstance(collision_tally, TallyCollision) assert tracklength_tally.spatial_filter_type == SPATIAL_FILTER_NONE @@ -28,3 +35,26 @@ def test_tally_factory_routing_invalid_score_mix(capsys): mcdc.Tally(scores=["flux", "energy_deposition"]) out = capsys.readouterr().out assert "Cannot mix tracklength scores with collision ones." in out + + +def test_tally_factory_routing_invalid_net_current_selectors(slab_plane_x, capsys): + with pytest.raises(SystemExit): + mcdc.Tally( + surface=slab_plane_x["s_mid"], + cell=slab_plane_x["c_right"], + scores=["net-current"], + ) + out = capsys.readouterr().out + assert "must specify exactly one of surface or cell" in out + + with pytest.raises(SystemExit): + mcdc.Tally(scores=["net-current"]) + out = capsys.readouterr().out + assert "needs either a surface or a cell tally" in out + + +def test_surface_tally_rejects_current_in_out_scores(slab_plane_x, capsys): + with pytest.raises(SystemExit): + mcdc.Tally(surface=slab_plane_x["s_mid"], scores=["current-in"]) + out = capsys.readouterr().out + assert "Surface tally currently supports only scores=['net-current']." in out From 48513310e6d8a205cc7636475b59e3bb277f1223 Mon Sep 17 00:00:00 2001 From: Massimo Larsen Date: Wed, 20 May 2026 09:26:34 -0700 Subject: [PATCH 03/15] Added better testing with more edge cases. Optimized cell tally. Fixed random issues --- mcdc/constant.py | 6 +- mcdc/object_/tally.py | 22 ++--- mcdc/transport/geometry/interface.py | 49 ++++++---- mcdc/transport/tally/score.py | 2 +- test/unit/transport/tally/conftest.py | 23 +++++ test/unit/transport/tally/test_cell_tally.py | 98 ++++++++++++------- .../tally/test_nonbounded_surface_tally.py | 43 ++++++++ .../tally/test_tally_factory_routing.py | 2 +- 8 files changed, 175 insertions(+), 70 deletions(-) diff --git a/mcdc/constant.py b/mcdc/constant.py index 27d31c8d3..9bd70766b 100644 --- a/mcdc/constant.py +++ b/mcdc/constant.py @@ -5,9 +5,9 @@ # Tally types TALLY_SURFACE = 0 -TALLY_CELL = 1 -TALLY_COLLISION = 2 -TALLY_TRACKLENGTH = 3 +TALLY_COLLISION = 1 +TALLY_TRACKLENGTH = 2 +TALLY_CELL = 3 # Tally spatial filter types SPATIAL_FILTER_NONE = 0 diff --git a/mcdc/object_/tally.py b/mcdc/object_/tally.py index 087b4e2b9..47826bcd2 100644 --- a/mcdc/object_/tally.py +++ b/mcdc/object_/tally.py @@ -50,10 +50,10 @@ from mcdc.object_.simulation import simulation from mcdc.print_ import print_1d_array, print_error -SURFACE_SCORES = set(["net-current"]) -CELL_CURRENT_SCORES = set(["net-current", "current-in", "current-out"]) -TRACKLENGTH_SCORES = set(["flux", "density", "collision", "capture", "fission"]) -COLLISION_SCORES = set(["energy_deposition"]) +SURFACE_SCORES = {"net-current"} +CELL_CURRENT_SCORES = {"net-current", "current-in", "current-out"} +TRACKLENGTH_SCORES = {"flux", "density", "collision", "capture", "fission"} +COLLISION_SCORES = {"energy_deposition"} class Tally(ObjectPolymorphic): @@ -105,25 +105,25 @@ def __new__( # Determine type and create the tally self based on the provided # spatial filters and scores - has_cell_current_score = any(score in CELL_CURRENT_SCORES for score in scores) + has_current_score = any(score in CELL_CURRENT_SCORES for score in scores) - # Surface/cell net-current tally - if has_cell_current_score: + # Surface/cell current tally + if has_current_score: for score in scores: if score not in CELL_CURRENT_SCORES: print_error( - "Cannot mix cell-current scores with non-current scores " - f"in one tally. Cell-current scores: {CELL_CURRENT_SCORES}." + "Cannot mix current scores with non-current scores " + f"in one tally. Current scores: {CELL_CURRENT_SCORES}." ) if surface is not None and cell is not None: print_error( - "Net-current tally must specify exactly one of surface or cell." + "Current tally must specify exactly one of surface or cell." ) if surface is None and cell is None: print_error( - "Scoring 'net-current' needs either a surface or a cell tally." + "Current scores need either a surface or a cell tally." ) if surface is not None: diff --git a/mcdc/transport/geometry/interface.py b/mcdc/transport/geometry/interface.py index 6268966e0..4b9d049dc 100644 --- a/mcdc/transport/geometry/interface.py +++ b/mcdc/transport/geometry/interface.py @@ -443,13 +443,8 @@ def surface_crossing(P_arr, simulation, data): P = P_arr[0] crossed_surface_ID = P["surface_ID"] - # Determine pre/post top-cell IDs for non-reflective transfers. - pre_cell_ID = -1 - post_cell_ID = -1 surface = simulation["surfaces"][crossed_surface_ID] BC = surface["boundary_condition"] - if BC != BC_REFLECTIVE: - pre_cell_ID, post_cell_ID = _get_crossing_top_cell_IDs(P_arr, simulation, data) # Apply BC if BC == BC_VACUUM: @@ -472,14 +467,17 @@ def surface_crossing(P_arr, simulation, data): # Score cell net-current tallies tied to this crossed surface. if BC != BC_REFLECTIVE: - for i in range(simulation["N_cell_tally"]): - tally = simulation["cell_tallies"][i] - cell = simulation["cells"][tally["cell_ID"]] - if not _cell_has_surface(cell, crossed_surface_ID, data): - continue - tally_module.score.cell_tally( - P_arr, tally, pre_cell_ID, post_cell_ID, simulation, data + if simulation["N_cell_tally"] > 0: + pre_cell_ID, post_cell_ID = _get_crossing_top_cell_IDs( + P_arr, simulation, data + ) + _score_cell_current_tallies( + P_arr, pre_cell_ID, pre_cell_ID, post_cell_ID, simulation, data ) + if post_cell_ID != pre_cell_ID: + _score_cell_current_tallies( + P_arr, post_cell_ID, pre_cell_ID, post_cell_ID, simulation, data + ) # Need to check new cell later? if P["alive"] and not BC == BC_REFLECTIVE: @@ -517,11 +515,28 @@ def _get_crossing_top_cell_IDs(particle_container, simulation, data): @njit -def _cell_has_surface(cell, surface_ID, data): - for i in range(cell["N_surface"]): - if int(mcdc_get.cell.surface_IDs(i, cell, data)) == surface_ID: - return True - return False +def _score_cell_current_tallies( + particle_container, + cell_ID, + pre_cell_ID, + post_cell_ID, + simulation, + data, +): + if cell_ID < 0: + return + + cell = simulation["cells"][cell_ID] + for i in range(cell["N_tally"]): + tally_base_ID = int(mcdc_get.cell.tally_IDs(i, cell, data)) + tally_base = simulation["tallies"][tally_base_ID] + if tally_base["child_type"] != TALLY_CELL: + continue + + tally = simulation["cell_tallies"][tally_base["child_ID"]] + tally_module.score.cell_tally( + particle_container, tally, pre_cell_ID, post_cell_ID, simulation, data + ) @njit diff --git a/mcdc/transport/tally/score.py b/mcdc/transport/tally/score.py index 8a08deca4..ccf6e305e 100644 --- a/mcdc/transport/tally/score.py +++ b/mcdc/transport/tally/score.py @@ -129,7 +129,7 @@ def cell_tally( score = particle["w"] elif score_type == SCORE_CURRENT_OUT: if exited: - score = -particle["w"] + score = particle["w"] if score != 0.0: util.atomic_add(data, idx_base + i_score, score) diff --git a/test/unit/transport/tally/conftest.py b/test/unit/transport/tally/conftest.py index 317676b85..ca58fd7e9 100644 --- a/test/unit/transport/tally/conftest.py +++ b/test/unit/transport/tally/conftest.py @@ -27,6 +27,29 @@ def material_mg(): return mcdc.MaterialMG(capture=np.array([1.0])) +@pytest.fixture +def crossing_particle(): + def _particle(surface_ID, x, ux, w=2.0): + particle_container = np.zeros(1, type_.particle) + particle = particle_container[0] + particle["alive"] = True + particle["particle_type"] = PARTICLE_NEUTRON + particle["surface_ID"] = surface_ID + particle["material_ID"] = 0 + particle["g"] = 0 + particle["x"] = x + particle["y"] = 0.0 + particle["z"] = 0.0 + particle["t"] = 0.0 + particle["ux"] = ux + particle["uy"] = 0.0 + particle["uz"] = 0.0 + particle["w"] = w + return particle_container + + return _particle + + @pytest.fixture def slab_plane_x(material_mg): # Two-cell 1D slab with shared interior PlaneX. diff --git a/test/unit/transport/tally/test_cell_tally.py b/test/unit/transport/tally/test_cell_tally.py index 2726ad383..e4aeb80d5 100644 --- a/test/unit/transport/tally/test_cell_tally.py +++ b/test/unit/transport/tally/test_cell_tally.py @@ -1,31 +1,10 @@ import numpy as np import mcdc -import mcdc.numba_types as type_ -from mcdc.constant import PARTICLE_NEUTRON from mcdc.main import preparation from mcdc.transport.geometry import interface as geometry_interface -def _new_particle(surface_ID, x, ux, w=2.0): - particle_container = np.zeros(1, type_.particle) - particle = particle_container[0] - particle["alive"] = True - particle["particle_type"] = PARTICLE_NEUTRON - particle["surface_ID"] = surface_ID - particle["material_ID"] = 0 - particle["g"] = 0 - particle["x"] = x - particle["y"] = 0.0 - particle["z"] = 0.0 - particle["t"] = 0.0 - particle["ux"] = ux - particle["uy"] = 0.0 - particle["uz"] = 0.0 - particle["w"] = w - return particle_container - - def _bin_value(cell_tally, mcdc_struct, data): tally_base = mcdc_struct["tallies"][cell_tally["parent_ID"]] return data[tally_base["bin_offset"]] @@ -36,24 +15,24 @@ def _bin_value_score(cell_tally, mcdc_struct, data, score_idx): return data[tally_base["bin_offset"] + score_idx] -def test_cell_tally_sign_for_incoming_and_outgoing(slab_plane_x): +def test_cell_tally_sign_for_incoming_and_outgoing(slab_plane_x, crossing_particle): cell_tally_obj = mcdc.Tally(cell=slab_plane_x["c_right"], scores=["net-current"]) mcdc_container, data = preparation() mcdc_struct = mcdc_container[0] cell_tally = mcdc_struct["cell_tallies"][cell_tally_obj.child_ID] # Left -> right across the shared interior surface: incoming to c_right (+) - particle_container = _new_particle(slab_plane_x["s_mid"].ID, x=0.0, ux=0.5) + particle_container = crossing_particle(slab_plane_x["s_mid"].ID, x=0.0, ux=0.5) geometry_interface.surface_crossing(particle_container, mcdc_struct, data) assert np.isclose(_bin_value(cell_tally, mcdc_struct, data), 2.0) # Right -> left across the shared interior surface: outgoing from c_right (-) - particle_container = _new_particle(slab_plane_x["s_mid"].ID, x=0.0, ux=-0.5) + particle_container = crossing_particle(slab_plane_x["s_mid"].ID, x=0.0, ux=-0.5) geometry_interface.surface_crossing(particle_container, mcdc_struct, data) assert np.isclose(_bin_value(cell_tally, mcdc_struct, data), 0.0) -def test_cell_tally_records_in_and_out_separately(slab_plane_x): +def test_cell_tally_records_in_and_out_separately(slab_plane_x, crossing_particle): cell_tally_obj = mcdc.Tally( cell=slab_plane_x["c_right"], scores=["net-current", "current-in", "current-out"], @@ -63,40 +42,85 @@ def test_cell_tally_records_in_and_out_separately(slab_plane_x): cell_tally = mcdc_struct["cell_tallies"][cell_tally_obj.child_ID] # One incoming and one outgoing crossing. - particle_container = _new_particle(slab_plane_x["s_mid"].ID, x=0.0, ux=0.5) + particle_container = crossing_particle(slab_plane_x["s_mid"].ID, x=0.0, ux=0.5) geometry_interface.surface_crossing(particle_container, mcdc_struct, data) - particle_container = _new_particle(slab_plane_x["s_mid"].ID, x=0.0, ux=-0.5) + particle_container = crossing_particle(slab_plane_x["s_mid"].ID, x=0.0, ux=-0.5) geometry_interface.surface_crossing(particle_container, mcdc_struct, data) # Scores are in requested order: [net-current, current-in, current-out]. - # current-out uses negative sign for exiting particles. + # Partial currents are positive; net-current carries the sign. assert np.isclose(_bin_value_score(cell_tally, mcdc_struct, data, 0), 0.0) assert np.isclose(_bin_value_score(cell_tally, mcdc_struct, data, 1), 2.0) - assert np.isclose(_bin_value_score(cell_tally, mcdc_struct, data, 2), -2.0) + assert np.isclose(_bin_value_score(cell_tally, mcdc_struct, data, 2), 2.0) -def test_cell_tally_counts_outgoing_to_vacuum(slab_plane_x): - cell_tally_obj = mcdc.Tally(cell=slab_plane_x["c_right"], scores=["net-current"]) +def test_cell_tally_counts_outgoing_to_vacuum(slab_plane_x, crossing_particle): + cell_tally_obj = mcdc.Tally( + cell=slab_plane_x["c_right"], + scores=["net-current", "current-in", "current-out"], + ) mcdc_container, data = preparation() mcdc_struct = mcdc_container[0] cell_tally = mcdc_struct["cell_tallies"][cell_tally_obj.child_ID] # c_right -> vacuum across the right boundary: outgoing (-) - particle_container = _new_particle(slab_plane_x["s_right"].ID, x=1.0, ux=0.5) + particle_container = crossing_particle(slab_plane_x["s_right"].ID, x=1.0, ux=0.5) + particle = particle_container[0] geometry_interface.surface_crossing(particle_container, mcdc_struct, data) - assert np.isclose(_bin_value(cell_tally, mcdc_struct, data), -2.0) + assert not particle["alive"] + assert np.isclose(_bin_value_score(cell_tally, mcdc_struct, data, 0), -2.0) + assert np.isclose(_bin_value_score(cell_tally, mcdc_struct, data, 1), 0.0) + assert np.isclose(_bin_value_score(cell_tally, mcdc_struct, data, 2), 2.0) -def test_cell_tally_ignores_reflective_crossing(material_mg): +def test_cell_tally_ignores_reflective_crossing(material_mg, crossing_particle): s_left = mcdc.Surface.PlaneX(x=-1.0, boundary_condition="vacuum") s_right = mcdc.Surface.PlaneX(x=1.0, boundary_condition="reflective") c_mid = mcdc.Cell(region=+s_left & -s_right, fill=material_mg) - cell_tally_obj = mcdc.Tally(cell=c_mid, scores=["net-current"]) + cell_tally_obj = mcdc.Tally( + cell=c_mid, + scores=["net-current", "current-in", "current-out"], + ) mcdc_container, data = preparation() mcdc_struct = mcdc_container[0] cell_tally = mcdc_struct["cell_tallies"][cell_tally_obj.child_ID] - particle_container = _new_particle(s_right.ID, x=1.0, ux=0.5) + particle_container = crossing_particle(s_right.ID, x=1.0, ux=0.5) + particle = particle_container[0] geometry_interface.surface_crossing(particle_container, mcdc_struct, data) - assert np.isclose(_bin_value(cell_tally, mcdc_struct, data), 0.0) + assert particle["alive"] + assert np.isclose(particle["ux"], -0.5) + assert np.isclose(_bin_value_score(cell_tally, mcdc_struct, data, 0), 0.0) + assert np.isclose(_bin_value_score(cell_tally, mcdc_struct, data, 1), 0.0) + assert np.isclose(_bin_value_score(cell_tally, mcdc_struct, data, 2), 0.0) + + +def test_cell_tally_scores_curved_boundary(material_mg, crossing_particle): + s_cyl = mcdc.Surface.CylinderZ(center=(0.0, 0.0), radius=1.0) + c_inner = mcdc.Cell(region=-s_cyl, fill=material_mg) + mcdc.Cell(region=+s_cyl, fill=material_mg) + cell_tally_obj = mcdc.Tally( + cell=c_inner, + scores=["net-current", "current-in", "current-out"], + ) + + mcdc_container, data = preparation() + mcdc_struct = mcdc_container[0] + cell_tally = mcdc_struct["cell_tallies"][cell_tally_obj.child_ID] + + # Inner -> outer across the cylindrical surface: outgoing from c_inner (-) + particle_container = crossing_particle(s_cyl.ID, x=1.0, ux=0.5) + geometry_interface.surface_crossing(particle_container, mcdc_struct, data) + + assert np.isclose(_bin_value_score(cell_tally, mcdc_struct, data, 0), -2.0) + assert np.isclose(_bin_value_score(cell_tally, mcdc_struct, data, 1), 0.0) + assert np.isclose(_bin_value_score(cell_tally, mcdc_struct, data, 2), 2.0) + + # Outer -> inner across the same curved surface: incoming to c_inner (+) + particle_container = crossing_particle(s_cyl.ID, x=1.0, ux=-0.5) + geometry_interface.surface_crossing(particle_container, mcdc_struct, data) + + assert np.isclose(_bin_value_score(cell_tally, mcdc_struct, data, 0), 0.0) + assert np.isclose(_bin_value_score(cell_tally, mcdc_struct, data, 1), 2.0) + assert np.isclose(_bin_value_score(cell_tally, mcdc_struct, data, 2), 2.0) diff --git a/test/unit/transport/tally/test_nonbounded_surface_tally.py b/test/unit/transport/tally/test_nonbounded_surface_tally.py index c6ba9c855..a4b4fc00e 100644 --- a/test/unit/transport/tally/test_nonbounded_surface_tally.py +++ b/test/unit/transport/tally/test_nonbounded_surface_tally.py @@ -1,6 +1,8 @@ import numpy as np import pytest +import mcdc +from mcdc.main import preparation from mcdc.transport.geometry import interface as geometry_interface @@ -52,3 +54,44 @@ def test_unbounded_surface_tally_scoring( assert np.isclose(bin_value(unbounded_tally, mcdc_struct, data), expected_unbounded) assert np.isclose(bin_value(bounded_tally, mcdc_struct, data), expected_bounded) + + +def test_surface_tally_scores_vacuum_boundary(material_mg, bin_value, crossing_particle): + s_left = mcdc.Surface.PlaneX(x=-1.0, boundary_condition="vacuum") + s_right = mcdc.Surface.PlaneX(x=1.0, boundary_condition="vacuum") + mcdc.Cell(region=+s_left & -s_right, fill=material_mg) + tally_obj = mcdc.Tally(surface=s_right, scores=["net-current"]) + + mcdc_container, data = preparation() + mcdc_struct = mcdc_container[0] + tally = mcdc_struct["surface_tallies"][tally_obj.child_ID] + + particle_container = crossing_particle(s_right.ID, x=1.0, ux=0.5) + particle = particle_container[0] + + geometry_interface.surface_crossing(particle_container, mcdc_struct, data) + + assert not particle["alive"] + assert np.isclose(bin_value(tally, mcdc_struct, data), 2.0) + + +def test_surface_tally_scores_after_reflective_boundary( + material_mg, bin_value, crossing_particle +): + s_left = mcdc.Surface.PlaneX(x=-1.0, boundary_condition="vacuum") + s_right = mcdc.Surface.PlaneX(x=1.0, boundary_condition="reflective") + mcdc.Cell(region=+s_left & -s_right, fill=material_mg) + tally_obj = mcdc.Tally(surface=s_right, scores=["net-current"]) + + mcdc_container, data = preparation() + mcdc_struct = mcdc_container[0] + tally = mcdc_struct["surface_tallies"][tally_obj.child_ID] + + particle_container = crossing_particle(s_right.ID, x=1.0, ux=0.5) + particle = particle_container[0] + + geometry_interface.surface_crossing(particle_container, mcdc_struct, data) + + assert particle["alive"] + assert np.isclose(particle["ux"], -0.5) + assert np.isclose(bin_value(tally, mcdc_struct, data), -2.0) diff --git a/test/unit/transport/tally/test_tally_factory_routing.py b/test/unit/transport/tally/test_tally_factory_routing.py index d033b1580..d922d324d 100644 --- a/test/unit/transport/tally/test_tally_factory_routing.py +++ b/test/unit/transport/tally/test_tally_factory_routing.py @@ -50,7 +50,7 @@ def test_tally_factory_routing_invalid_net_current_selectors(slab_plane_x, capsy with pytest.raises(SystemExit): mcdc.Tally(scores=["net-current"]) out = capsys.readouterr().out - assert "needs either a surface or a cell tally" in out + assert "Current scores need either a surface or a cell tally" in out def test_surface_tally_rejects_current_in_out_scores(slab_plane_x, capsys): From cdd74c973ae047c38a36be1ef990fbbe90796808 Mon Sep 17 00:00:00 2001 From: Massimo Larsen Date: Wed, 20 May 2026 09:49:33 -0700 Subject: [PATCH 04/15] ran black formatting --- mcdc/code_factory/numba_objects_generator.py | 6 ++++-- mcdc/object_/tally.py | 4 +--- test/unit/transport/tally/test_nonbounded_surface_tally.py | 4 +++- 3 files changed, 8 insertions(+), 6 deletions(-) diff --git a/mcdc/code_factory/numba_objects_generator.py b/mcdc/code_factory/numba_objects_generator.py index c466b92d6..1ee8d7306 100644 --- a/mcdc/code_factory/numba_objects_generator.py +++ b/mcdc/code_factory/numba_objects_generator.py @@ -831,8 +831,10 @@ def align(field_list): pad_id = 0 for field in field_list: if len(field) > 3: - print_error("Unexpected struct field specification. Specifications \ - usually only consist of 3 or fewer members") + print_error( + "Unexpected struct field specification. Specifications \ + usually only consist of 3 or fewer members" + ) multiplier = 1 if len(field) == 3: field = (field[0], field[1], fixup_dims(field[2])) diff --git a/mcdc/object_/tally.py b/mcdc/object_/tally.py index 47826bcd2..f4602509c 100644 --- a/mcdc/object_/tally.py +++ b/mcdc/object_/tally.py @@ -122,9 +122,7 @@ def __new__( ) if surface is None and cell is None: - print_error( - "Current scores need either a surface or a cell tally." - ) + print_error("Current scores need either a surface or a cell tally.") if surface is not None: if not set(scores) <= SURFACE_SCORES: diff --git a/test/unit/transport/tally/test_nonbounded_surface_tally.py b/test/unit/transport/tally/test_nonbounded_surface_tally.py index a4b4fc00e..03e59c51e 100644 --- a/test/unit/transport/tally/test_nonbounded_surface_tally.py +++ b/test/unit/transport/tally/test_nonbounded_surface_tally.py @@ -56,7 +56,9 @@ def test_unbounded_surface_tally_scoring( assert np.isclose(bin_value(bounded_tally, mcdc_struct, data), expected_bounded) -def test_surface_tally_scores_vacuum_boundary(material_mg, bin_value, crossing_particle): +def test_surface_tally_scores_vacuum_boundary( + material_mg, bin_value, crossing_particle +): s_left = mcdc.Surface.PlaneX(x=-1.0, boundary_condition="vacuum") s_right = mcdc.Surface.PlaneX(x=1.0, boundary_condition="vacuum") mcdc.Cell(region=+s_left & -s_right, fill=material_mg) From f86e9938c4e31118722dc1c066ef86239c33f266 Mon Sep 17 00:00:00 2001 From: Massimo Larsen Date: Wed, 20 May 2026 13:24:52 -0700 Subject: [PATCH 05/15] Initial test update commit. Renamed tests, modified distribution tests to run with numba, modified surface tests to run with pytest. --- pyproject.toml | 1 + test/unit/conftest.py | 22 +++++ ...eposition.py => test_energy_deposition.py} | 0 test/unit/transport/distributions/conftest.py | 90 +++++++++++-------- .../distributions/test_evaporation.py | 23 +++-- .../distributions/test_kalbach_mann.py | 13 +-- .../distributions/test_level_scattering.py | 7 +- .../distributions/test_maxwellian.py | 23 +++-- .../distributions/test_multi_table.py | 23 +++-- .../distributions/test_nbody_correlated.py | 43 ++++++--- .../test_tabulated_distribution.py | 17 ++-- .../test_tabulated_energy_angle.py | 14 ++- .../transport/geometry/surface/conftest.py | 25 ++++++ .../surface/{plane_x.py => test_plane_x.py} | 44 ++++----- .../surface/{plane_y.py => test_plane_y.py} | 44 ++++----- .../surface/{plane_z.py => test_plane_z.py} | 44 ++++----- .../surface/{torus_z.py => test_torus_z.py} | 44 ++++----- test/unit/transport/tally/conftest.py | 5 -- ...ight_windows.py => test_weight_windows.py} | 0 .../util/{find_bin.py => test_find_bin.py} | 0 ...lation.py => test_linear_interpolation.py} | 0 21 files changed, 307 insertions(+), 175 deletions(-) create mode 100644 test/unit/conftest.py rename test/unit/object_/tally/{energy_deposition.py => test_energy_deposition.py} (100%) create mode 100644 test/unit/transport/geometry/surface/conftest.py rename test/unit/transport/geometry/surface/{plane_x.py => test_plane_x.py} (96%) rename test/unit/transport/geometry/surface/{plane_y.py => test_plane_y.py} (96%) rename test/unit/transport/geometry/surface/{plane_z.py => test_plane_z.py} (96%) rename test/unit/transport/geometry/surface/{torus_z.py => test_torus_z.py} (97%) rename test/unit/transport/technique/{weight_windows.py => test_weight_windows.py} (100%) rename test/unit/transport/util/{find_bin.py => test_find_bin.py} (100%) rename test/unit/transport/util/{linear_interpolation.py => test_linear_interpolation.py} (100%) diff --git a/pyproject.toml b/pyproject.toml index 6e7c63049..2a7c9850c 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -95,5 +95,6 @@ force-exclude = ''' mcdc/numba_types\.py | mcdc/mcdc_get/ | mcdc/mcdc_set/ + | mcdc/lib/ ) ''' diff --git a/test/unit/conftest.py b/test/unit/conftest.py new file mode 100644 index 000000000..6ae1f6130 --- /dev/null +++ b/test/unit/conftest.py @@ -0,0 +1,22 @@ +import os + + +def pytest_addoption(parser): + parser.addoption( + "--mode", + choices=["python", "numba"], + default="python", + help="MCDC execution mode for unit tests.", + ) + + +def pytest_configure(config): + mode = config.getoption("--mode") + if mode == "python": + os.environ["NUMBA_DISABLE_JIT"] = "1" + else: + os.environ.pop("NUMBA_DISABLE_JIT", None) + + +def pytest_report_header(config): + return f"MCDC mode: {config.getoption('--mode')}" diff --git a/test/unit/object_/tally/energy_deposition.py b/test/unit/object_/tally/test_energy_deposition.py similarity index 100% rename from test/unit/object_/tally/energy_deposition.py rename to test/unit/object_/tally/test_energy_deposition.py diff --git a/test/unit/transport/distributions/conftest.py b/test/unit/transport/distributions/conftest.py index 96ee1b6c9..994845027 100644 --- a/test/unit/transport/distributions/conftest.py +++ b/test/unit/transport/distributions/conftest.py @@ -1,54 +1,66 @@ -import os +import numpy as np import pytest - -# Force pure-Python execution for numba in tests. -os.environ.setdefault("NUMBA_DISABLE_JIT", "1") +from numba import njit import mcdc.transport.distribution as dist +MOCK_RNG_MAX_VALUES = 8 +MOCK_RNG_STATE_DTYPE = np.dtype( + [ + ("idx", np.int64), + ("n_values", np.int64), + ("values", np.float64, (MOCK_RNG_MAX_VALUES,)), + ] +) -class MockRNG: - def __init__(self, values): - self._values = list(values) - self._i = 0 - def lcg(self, _state_container): - # MCDC uses dist.rng.lcg(state) as its RNG hook. We override it with this - # deterministic sequence so tests can be analytic and reproducible. - if self._i >= len(self._values): - raise IndexError("MockRNG depleted") - value = self._values[self._i] - self._i += 1 - return value +@njit +def mock_lcg_known_sequence(rng_state): + # Force lcg to pick mock rng values in sequence + state = rng_state[0] + i = state["idx"] + n = state["n_values"] + if i < n: + value = state["values"][i] + else: + value = np.nan -@pytest.fixture -def rng_state(): - """ - Return the minimal RNG state container expected by MCDC. + state["idx"] = i + 1 + return value + + +def mock_rng(*values): + if len(values) == 1 and isinstance(values[0], (list, tuple, np.ndarray)): + values = tuple(values[0]) - MCDC passes a list of per-thread state dicts; tests only need one. - """ - return [dict(rng_seed=0)] + state = np.zeros(1, dtype=MOCK_RNG_STATE_DTYPE) + state[0]["idx"] = 0 + state[0]["n_values"] = len(values) + if values: + state[0]["values"][: len(values)] = np.asarray(values, dtype=np.float64) + return state @pytest.fixture -def rng_sequence(): - """ - Temporarily replace dist.rng.lcg with a deterministic sequence. - - Usage: - rng_sequence([0.1, 0.2, 0.3]) - ... code under test that calls dist.rng.lcg(...) - """ - original_lcg = dist.rng.lcg +def make_distribution_record(): + # Build one typed record from a dict of field values. + def _make(record_dtype, values_dict): + container = np.zeros(1, record_dtype) + record = container[0] + for key, value in values_dict.items(): + record[key] = value + return record + + return _make - def _apply(values): - # Swap in a mock LCG implementation that returns the provided values. - rng = MockRNG(values) - dist.rng.lcg = rng.lcg - return rng - # Yield a callable to the test, then restore the real RNG after the test ends. - yield _apply +@pytest.fixture +def mock_rng_sequence(): + # Replace mcdc lcg with the known sequence + original_lcg = dist.rng.lcg + dist.rng.lcg = mock_lcg_known_sequence + # Yield the mock_rng sequence + yield mock_rng + # Reset the lcg. If the test fails, the mock_rng could still be installed and cause later tests to fail dist.rng.lcg = original_lcg diff --git a/test/unit/transport/distributions/test_evaporation.py b/test/unit/transport/distributions/test_evaporation.py index 354d131b6..9a5774133 100644 --- a/test/unit/transport/distributions/test_evaporation.py +++ b/test/unit/transport/distributions/test_evaporation.py @@ -1,20 +1,31 @@ import math +import numpy as np +import mcdc.numba_types as type_ import mcdc.transport.distribution as dist from .test_data import make_test_table_data_constant -def test_evaporation_sample(rng_sequence, rng_state): +def test_evaporation_sample(mock_rng_sequence, make_distribution_record): # MCNP Theory & User Manual §2.4.3.5.4.7 (Law 9: Evaporation Spectrum) - table, data = make_test_table_data_constant(1.0) - mcdc = {"table_data": [table]} - evaporation = {"nuclear_temperature_ID": 0, "restriction_energy": 0.0} + table_dict, data = make_test_table_data_constant(1.0) + table = make_distribution_record(type_.table_data, table_dict) + evaporation = make_distribution_record( + type_.evaporation_distribution, + {"nuclear_temperature_ID": 0, "restriction_energy": 0.0}, + ) + data = np.asarray(data, dtype=np.float64) + + simulation_dtype = np.dtype([("table_data", type_.table_data, (1,))]) + simulation_container = np.zeros(1, dtype=simulation_dtype) + simulation = simulation_container[0] + simulation["table_data"][0] = table xi1, xi2 = 0.1, 0.2 - rng_sequence([xi1, xi2]) + mock_rng = mock_rng_sequence(xi1, xi2) - sampled_E = dist.sample_evaporation(2.0, rng_state, evaporation, mcdc, data) + sampled_E = dist.sample_evaporation(2.0, mock_rng, evaporation, simulation, data) # The constant temperature table gives T(E_in) = 1.0 for this test. # With U = 0, the MCNP notation E_in - U becomes 2.0, so # w = (E_in - U) / T(E_in) = 2.0. diff --git a/test/unit/transport/distributions/test_kalbach_mann.py b/test/unit/transport/distributions/test_kalbach_mann.py index 2d29048c3..51f912b6f 100644 --- a/test/unit/transport/distributions/test_kalbach_mann.py +++ b/test/unit/transport/distributions/test_kalbach_mann.py @@ -1,20 +1,23 @@ import math +import numpy as np +import mcdc.numba_types as type_ import mcdc.transport.distribution as dist from .test_data import make_test_kalbach_mann_data -def test_kalbach_mann_sample(rng_sequence, rng_state): +def test_kalbach_mann_sample(mock_rng_sequence, make_distribution_record): # MCNP Theory & User Manual §2.4.3.5.4.11 (Law 44: Kalbach-87 Correlated Energy-angle Scattering) - kalbach, data = make_test_kalbach_mann_data() + kalbach_dict, data = make_test_kalbach_mann_data() + kalbach = make_distribution_record(type_.kalbach_mann_distribution, kalbach_dict) + data = np.asarray(data, dtype=np.float64) # For E_in = 2.0 on the grid [1, 3], the interpolation fraction is r = 0.5. # xi_1 = 0.3 < r, so the second energy table is selected. xi1, xi2, xi3, xi4 = 0.3, 0.1, 0.7, 0.5 - rng_sequence([xi1, xi2, xi3, xi4]) - - sampled_E, sampled_mu = dist.sample_kalbach_mann(2.0, rng_state, kalbach, data) + mock_rng = mock_rng_sequence(xi1, xi2, xi3, xi4) + sampled_E, sampled_mu = dist.sample_kalbach_mann(2.0, mock_rng, kalbach, data) # Law 44 uses the Law 4 energy construction, so with xi_2 = 0.1 in the first bin # of the selected table: diff --git a/test/unit/transport/distributions/test_level_scattering.py b/test/unit/transport/distributions/test_level_scattering.py index 2ecba06e6..45ca3bf6a 100644 --- a/test/unit/transport/distributions/test_level_scattering.py +++ b/test/unit/transport/distributions/test_level_scattering.py @@ -1,15 +1,18 @@ import math +import mcdc.numba_types as type_ import mcdc.transport.distribution as dist -def test_level_scattering_sample(): +def test_level_scattering_sample(make_distribution_record): # MCNP Theory & User Manual §2.4.3.5.4.3 (Law 3: Inelastic Scattering from Nuclear Levels) # The implementation stores the Law 3 relation in the linear form # E_out = C2 * (E_in - C1). # For this test, E_in = 5, C1 = 1, and C2 = 0.5, so # E_out = 0.5 * (5 - 1) = 2. - level = {"C1": 1.0, "C2": 0.5} + level = make_distribution_record( + type_.level_scattering_distribution, {"C1": 1.0, "C2": 0.5} + ) sampled_E = dist.sample_level_scattering(5.0, level) expected_E = 2.0 assert math.isclose(sampled_E, expected_E, rel_tol=0.0, abs_tol=1e-12) diff --git a/test/unit/transport/distributions/test_maxwellian.py b/test/unit/transport/distributions/test_maxwellian.py index 74fa5bec2..44c85a354 100644 --- a/test/unit/transport/distributions/test_maxwellian.py +++ b/test/unit/transport/distributions/test_maxwellian.py @@ -1,21 +1,32 @@ import math +import numpy as np +import mcdc.numba_types as type_ import mcdc.transport.distribution as dist from mcdc.constant import PI from .test_data import make_test_table_data_constant -def test_maxwellian_sample(rng_sequence, rng_state): +def test_maxwellian_sample(mock_rng_sequence, make_distribution_record): # MCNP Theory & User Manual §2.4.3.5.4.6 (Law 7: Simple Maxwell Fission Spectrum) - table, data = make_test_table_data_constant(1.0) - mcdc = {"table_data": [table]} - maxwellian = {"nuclear_temperature_ID": 0, "restriction_energy": 0.0} + table_dict, data = make_test_table_data_constant(1.0) + table = make_distribution_record(type_.table_data, table_dict) + maxwellian = make_distribution_record( + type_.maxwellian_distribution, + {"nuclear_temperature_ID": 0, "restriction_energy": 0.0}, + ) + data = np.asarray(data, dtype=np.float64) + + simulation_dtype = np.dtype([("table_data", type_.table_data, (1,))]) + simulation_container = np.zeros(1, dtype=simulation_dtype) + simulation = simulation_container[0] + simulation["table_data"][0] = table xi1, xi2, xi3 = 0.9, 0.9, 0.0 - rng_sequence([xi1, xi2, xi3]) + mock_rng = mock_rng_sequence(xi1, xi2, xi3) - sampled_E = dist.sample_maxwellian(2.0, rng_state, maxwellian, mcdc, data) + sampled_E = dist.sample_maxwellian(2.0, mock_rng, maxwellian, simulation, data) # MCNP Eq. (2.73): # E_out = -T(E_in) * [xi_1^2 / (xi_1^2 + xi_2^2) * ln(xi_3) + ln(xi_4)] # MCDC uses the equivalent polar-form reduction diff --git a/test/unit/transport/distributions/test_multi_table.py b/test/unit/transport/distributions/test_multi_table.py index d1a750758..bf64c65c1 100644 --- a/test/unit/transport/distributions/test_multi_table.py +++ b/test/unit/transport/distributions/test_multi_table.py @@ -1,24 +1,31 @@ import math +import numpy as np +import mcdc.numba_types as type_ import mcdc.transport.distribution as dist from .test_data import make_test_multi_table_data -def test_multi_table_distribution_sample(rng_sequence, rng_state): +def test_multi_table_distribution_sample(mock_rng_sequence, make_distribution_record): # MCNP Theory & User Manual §2.4.3.5.4.4 (Law 4: Tabular Distribution) - multi_table, data = make_test_multi_table_data() + multi_table_dict, data = make_test_multi_table_data() + multi_table = make_distribution_record( + type_.multi_table_distribution, multi_table_dict + ) + data = np.asarray(data, dtype=np.float64) # For E_in = 2.0 on the grid [1, 3], Eq. (2.62) gives r = 0.5. - # xi_2 = 0.3 < r, so Eq. (2.64) selects l = i + 1, i.e. the second table. - rng_sequence([0.3, 0.2]) + # xi_1 = 0.3 < r, so Eq. (2.64) selects l = i + 1, i.e. the second table. + xi1, xi2 = 0.3, 0.2 + mock_rng = mock_rng_sequence(xi1, xi2) - sampled_E = dist._sample_multi_table(2.0, rng_state, multi_table, data, scale=True) + sampled_E = dist._sample_multi_table(2.0, mock_rng, multi_table, data, scale=True) - # In the selected table, xi_1 = 0.2 falls in the first continuous bin. - # Eq. (2.65) gives E' = E_l,k + (xi_1 - c_l,k) / p_l,k = 100 + 0.2 / 0.01 = 120. + # In the selected table, xi_2 = 0.2 falls in the first continuous bin. + # Eq. (2.65) gives E' = E_l,k + (xi_2 - c_l,k) / p_l,k = 100 + 0.2 / 0.01 = 120. # The test data use constant p within the bin, so the linear-linear form in # Eq. (2.66) reduces to the same result. - E_prime = 100.0 + (0.2 - 0.0) / 0.01 + E_prime = 100.0 + (xi2 - 0.0) / 0.01 # Eq. (2.67) and Eq. (2.68) give the scaled bounds: # E_1 = 10 + 0.5 * (100 - 10) = 55 # E_K = 30 + 0.5 * (300 - 30) = 165 diff --git a/test/unit/transport/distributions/test_nbody_correlated.py b/test/unit/transport/distributions/test_nbody_correlated.py index a35f63361..dc2536783 100644 --- a/test/unit/transport/distributions/test_nbody_correlated.py +++ b/test/unit/transport/distributions/test_nbody_correlated.py @@ -1,25 +1,48 @@ import math +import numpy as np +import mcdc.numba_types as type_ import mcdc.transport.distribution as dist from mcdc.constant import DISTRIBUTION_N_BODY from .test_data import make_test_tabulated_data -def test_nbody_sample_correlated(rng_sequence, rng_state): +def test_nbody_sample_correlated(mock_rng_sequence, make_distribution_record): # MCNP Theory & User Manual §2.4.3.5.4.13 (Law 66: N-body Phase Space Distribution) - table, data = make_test_tabulated_data([2.0, 4.0, 6.0], [0.0, 0.4, 1.0]) - mcdc = {"nbody_distributions": [table]} - distribution = {"child_type": DISTRIBUTION_N_BODY, "child_ID": 0} + table_dict, data = make_test_tabulated_data([2.0, 4.0, 6.0], [0.0, 0.4, 1.0]) + nbody = make_distribution_record(type_.nbody_distribution, table_dict) + distribution = make_distribution_record( + type_.distribution, {"child_type": DISTRIBUTION_N_BODY, "child_ID": 0} + ) + data = np.asarray(data, dtype=np.float64) + + # Numba compiles all correlated-branch field accesses, so this container needs + # the three correlated distribution arrays even though this test uses N-body only. + simulation_dtype = np.dtype( + [ + ("kalbach_mann_distributions", type_.kalbach_mann_distribution, (1,)), + ( + "tabulated_energy_angle_distributions", + type_.tabulated_energy_angle_distribution, + (1,), + ), + ("nbody_distributions", type_.nbody_distribution, (1,)), + ] + ) + simulation_container = np.zeros(1, dtype=simulation_dtype) + simulation = simulation_container[0] + simulation["nbody_distributions"][0] = nbody # First value samples energy, second value samples isotropic cosine. - rng_sequence([0.2, 0.75]) + xi1, xi2 = 0.2, 0.75 + mock_rng = mock_rng_sequence(xi1, xi2) sampled_E, sampled_mu = dist.sample_correlated_distribution( 2.0, distribution, - rng_state, - mcdc, + mock_rng, + simulation, data, ) @@ -27,11 +50,11 @@ def test_nbody_sample_correlated(rng_sequence, rng_state): # samples the cosine isotropically. This test is therefore checking the current # reduced implementation, not reconstructing the full Law 66 rejection sampler # from Eq. (2.103) through Eq. (2.106). - # For the tabulated-energy part, xi = 0.2 gives linear interpolation in the first + # For the tabulated-energy part, xi_1 = 0.2 gives linear interpolation in the first # bin. For the angular part, MCNP Eq. (2.107) gives mu = 2 * xi_10 - 1 for # isotropic center-of-mass sampling. - expected_E = 2.0 + (0.2 - 0.0) * (4.0 - 2.0) / (0.4 - 0.0) - expected_mu = 2.0 * 0.75 - 1.0 + expected_E = 2.0 + (xi1 - 0.0) * (4.0 - 2.0) / (0.4 - 0.0) + expected_mu = 2.0 * xi2 - 1.0 assert math.isclose(sampled_E, expected_E, rel_tol=0.0, abs_tol=1e-12) assert math.isclose(sampled_mu, expected_mu, rel_tol=0.0, abs_tol=1e-12) diff --git a/test/unit/transport/distributions/test_tabulated_distribution.py b/test/unit/transport/distributions/test_tabulated_distribution.py index 648e1da4f..93e1d1166 100644 --- a/test/unit/transport/distributions/test_tabulated_distribution.py +++ b/test/unit/transport/distributions/test_tabulated_distribution.py @@ -1,19 +1,24 @@ import math +import numpy as np +import mcdc.numba_types as type_ import mcdc.transport.distribution as dist from .test_data import make_test_tabulated_data -def test_tabulated_distribution_sample(rng_sequence, rng_state): +def test_tabulated_distribution_sample(mock_rng_sequence, make_distribution_record): # MCNP Theory & User Manual §2.4.3.5.4.4 (Law 4: Tabular Distribution) - table, data = make_test_tabulated_data([1.0, 3.0, 7.0], [0.0, 0.4, 1.0]) - rng_sequence([0.2]) + table_dict, data = make_test_tabulated_data([1.0, 3.0, 7.0], [0.0, 0.4, 1.0]) + table = make_distribution_record(type_.tabulated_distribution, table_dict) + data = np.asarray(data, dtype=np.float64) + xi1 = 0.2 + mock_rng = mock_rng_sequence(xi1) - sampled_E = dist.sample_tabulated(table, rng_state, data) + sampled_E = dist.sample_tabulated(table, mock_rng, data) # This is the single-table inverse-CDF interpolation used by the tabulated sampler: - # xi = 0.2 lies in the first bin, so linear interpolation between + # xi_1 = 0.2 lies in the first bin, so linear interpolation between # (c_0, E_0) = (0.0, 1.0) and (c_1, E_1) = (0.4, 3.0) gives the expected value. - expected_E = 1.0 + (0.2 - 0.0) * (3.0 - 1.0) / (0.4 - 0.0) + expected_E = 1.0 + (xi1 - 0.0) * (3.0 - 1.0) / (0.4 - 0.0) assert math.isclose(sampled_E, expected_E, rel_tol=0.0, abs_tol=1e-12) diff --git a/test/unit/transport/distributions/test_tabulated_energy_angle.py b/test/unit/transport/distributions/test_tabulated_energy_angle.py index 8c7d770c7..8ed4a5263 100644 --- a/test/unit/transport/distributions/test_tabulated_energy_angle.py +++ b/test/unit/transport/distributions/test_tabulated_energy_angle.py @@ -1,19 +1,25 @@ import math +import numpy as np +import mcdc.numba_types as type_ import mcdc.transport.distribution as dist from .test_data import make_test_tabulated_energy_angle_data -def test_tabulated_energy_angle_sample(rng_sequence, rng_state): +def test_tabulated_energy_angle_sample(mock_rng_sequence, make_distribution_record): # MCNP Theory & User Manual §2.4.3.5.4.12 (Law 61: Correlated Energy-angle Scattering) - table, data = make_test_tabulated_energy_angle_data() + table_dict, data = make_test_tabulated_energy_angle_data() + table = make_distribution_record( + type_.tabulated_energy_angle_distribution, table_dict + ) + data = np.asarray(data, dtype=np.float64) xi1, xi2, xi3 = 0.3, 0.1, 0.25 - rng_sequence([xi1, xi2, xi3]) + mock_rng = mock_rng_sequence(xi1, xi2, xi3) sampled_E, sampled_mu = dist.sample_tabulated_energy_angle( - 2.0, rng_state, table, data + 2.0, mock_rng, table, data ) # Law 61 uses the Law 4 energy construction first. diff --git a/test/unit/transport/geometry/surface/conftest.py b/test/unit/transport/geometry/surface/conftest.py new file mode 100644 index 000000000..b26431eda --- /dev/null +++ b/test/unit/transport/geometry/surface/conftest.py @@ -0,0 +1,25 @@ +import numpy as np +import pytest + + +@pytest.fixture(autouse=True) +def reset_simulation(): + from mcdc.object_.simulation import simulation + + simulation.__init__() + yield + simulation.__init__() + + +@pytest.fixture +def compile_surfaces(): + def _compile(static_surface_obj, moving_surface_obj): + from mcdc.main import preparation + + structure_container, data = preparation() + structure = structure_container[0] + static_surface = structure["surfaces"][static_surface_obj.ID] + moving_surface = structure["surfaces"][moving_surface_obj.ID] + return static_surface, moving_surface, data + + return _compile diff --git a/test/unit/transport/geometry/surface/plane_x.py b/test/unit/transport/geometry/surface/test_plane_x.py similarity index 96% rename from test/unit/transport/geometry/surface/plane_x.py rename to test/unit/transport/geometry/surface/test_plane_x.py index 20fc2b365..f52694bdf 100644 --- a/test/unit/transport/geometry/surface/plane_x.py +++ b/test/unit/transport/geometry/surface/test_plane_x.py @@ -1,5 +1,6 @@ import mcdc import numpy as np +import pytest #### @@ -7,7 +8,7 @@ COINCIDENCE_TOLERANCE, INF, ) -from mcdc.main import preparation +import mcdc.numba_types as type_ # ====================================================================================== # Setup @@ -19,26 +20,11 @@ velocities = np.zeros((3, 3)) velocities[:, 0] = np.array([-1.0, 2.0, -3.0]) -# Test object: static surface -static_surface = mcdc.Surface.PlaneX(x=X) - -# Test object: moving surface -moving_surface = mcdc.Surface.PlaneX(x=X) -moving_surface.move(velocities, durations) - -# Create the dummy simulation structure and data -structure_container, data = preparation() -structure = structure_container[0] - -# Get the "compiled" test objects -static_surface = structure["surfaces"][0] -moving_surface = structure["surfaces"][1] - -# Particle object for testing -import mcdc.numba_types as type_ - -particle_container = np.zeros(1, type_.particle_data) -particle = particle_container[0] +static_surface = None +moving_surface = None +data = None +particle_container = None +particle = None # Miscellanies TINY = COINCIDENCE_TOLERANCE * 0.8 # Tiny value within coincidence tolerance @@ -49,6 +35,22 @@ plane_x, ) + +@pytest.fixture(autouse=True) +def setup_geometry_case(compile_surfaces): + global static_surface, moving_surface, data, particle_container, particle + + static_surface_obj = mcdc.Surface.PlaneX(x=X) + moving_surface_obj = mcdc.Surface.PlaneX(x=X) + moving_surface_obj.move(velocities, durations) + + static_surface, moving_surface, data = compile_surfaces( + static_surface_obj, moving_surface_obj + ) + particle_container = np.zeros(1, type_.particle_data) + particle = particle_container[0] + + # ===================================================================================== # Plane-X core functions # ===================================================================================== diff --git a/test/unit/transport/geometry/surface/plane_y.py b/test/unit/transport/geometry/surface/test_plane_y.py similarity index 96% rename from test/unit/transport/geometry/surface/plane_y.py rename to test/unit/transport/geometry/surface/test_plane_y.py index f9e1467d3..7f85629b9 100644 --- a/test/unit/transport/geometry/surface/plane_y.py +++ b/test/unit/transport/geometry/surface/test_plane_y.py @@ -1,5 +1,6 @@ import mcdc import numpy as np +import pytest #### @@ -7,7 +8,7 @@ COINCIDENCE_TOLERANCE, INF, ) -from mcdc.main import preparation +import mcdc.numba_types as type_ # ====================================================================================== # Setup @@ -19,26 +20,11 @@ velocities = np.zeros((3, 3)) velocities[:, 1] = np.array([-1.0, 2.0, -3.0]) -# Test object: static surface -static_surface = mcdc.Surface.PlaneY(y=Y) - -# Test object: moving surface -moving_surface = mcdc.Surface.PlaneY(y=Y) -moving_surface.move(velocities, durations) - -# Create the dummy simulation structure and data -structure_container, data = preparation() -structure = structure_container[0] - -# Get the "compiled" test objects -static_surface = structure["surfaces"][0] -moving_surface = structure["surfaces"][1] - -# Particle object for testing -import mcdc.numba_types as type_ - -particle_container = np.zeros(1, type_.particle_data) -particle = particle_container[0] +static_surface = None +moving_surface = None +data = None +particle_container = None +particle = None # Miscellanies TINY = COINCIDENCE_TOLERANCE * 0.8 # Tiny value within coincidence tolerance @@ -49,6 +35,22 @@ plane_y, ) + +@pytest.fixture(autouse=True) +def setup_geometry_case(compile_surfaces): + global static_surface, moving_surface, data, particle_container, particle + + static_surface_obj = mcdc.Surface.PlaneY(y=Y) + moving_surface_obj = mcdc.Surface.PlaneY(y=Y) + moving_surface_obj.move(velocities, durations) + + static_surface, moving_surface, data = compile_surfaces( + static_surface_obj, moving_surface_obj + ) + particle_container = np.zeros(1, type_.particle_data) + particle = particle_container[0] + + # ===================================================================================== # Plane-Y core functions # ===================================================================================== diff --git a/test/unit/transport/geometry/surface/plane_z.py b/test/unit/transport/geometry/surface/test_plane_z.py similarity index 96% rename from test/unit/transport/geometry/surface/plane_z.py rename to test/unit/transport/geometry/surface/test_plane_z.py index c4f591ee6..1d68016b2 100644 --- a/test/unit/transport/geometry/surface/plane_z.py +++ b/test/unit/transport/geometry/surface/test_plane_z.py @@ -1,5 +1,6 @@ import mcdc import numpy as np +import pytest #### @@ -7,7 +8,7 @@ COINCIDENCE_TOLERANCE, INF, ) -from mcdc.main import preparation +import mcdc.numba_types as type_ # ====================================================================================== # Setup @@ -19,26 +20,11 @@ velocities = np.zeros((3, 3)) velocities[:, 2] = np.array([-1.0, 2.0, -3.0]) -# Test object: static surface -static_surface = mcdc.Surface.PlaneZ(z=Z) - -# Test object: moving surface -moving_surface = mcdc.Surface.PlaneZ(z=Z) -moving_surface.move(velocities, durations) - -# Create the dummy simulation structure and data -structure_container, data = preparation() -structure = structure_container[0] - -# Get the "compiled" test objects -static_surface = structure["surfaces"][0] -moving_surface = structure["surfaces"][1] - -# Particle object for testing -import mcdc.numba_types as type_ - -particle_container = np.zeros(1, type_.particle_data) -particle = particle_container[0] +static_surface = None +moving_surface = None +data = None +particle_container = None +particle = None # Miscellanies TINY = COINCIDENCE_TOLERANCE * 0.8 # Tiny value within coincidence tolerance @@ -49,6 +35,22 @@ plane_z, ) + +@pytest.fixture(autouse=True) +def setup_geometry_case(compile_surfaces): + global static_surface, moving_surface, data, particle_container, particle + + static_surface_obj = mcdc.Surface.PlaneZ(z=Z) + moving_surface_obj = mcdc.Surface.PlaneZ(z=Z) + moving_surface_obj.move(velocities, durations) + + static_surface, moving_surface, data = compile_surfaces( + static_surface_obj, moving_surface_obj + ) + particle_container = np.zeros(1, type_.particle_data) + particle = particle_container[0] + + # ===================================================================================== # Plane-Z core functions # ===================================================================================== diff --git a/test/unit/transport/geometry/surface/torus_z.py b/test/unit/transport/geometry/surface/test_torus_z.py similarity index 97% rename from test/unit/transport/geometry/surface/torus_z.py rename to test/unit/transport/geometry/surface/test_torus_z.py index c3f83468d..cbbeba236 100644 --- a/test/unit/transport/geometry/surface/torus_z.py +++ b/test/unit/transport/geometry/surface/test_torus_z.py @@ -1,6 +1,7 @@ import mcdc import math import numpy as np +import pytest #### @@ -8,7 +9,7 @@ COINCIDENCE_TOLERANCE, INF, ) -from mcdc.main import preparation +import mcdc.numba_types as type_ # ====================================================================================== # Setup @@ -24,26 +25,11 @@ velocities = np.zeros((3, 3)) velocities[:, 0] = np.array([-1.0, 2.0, -3.0]) -# Test object: static surface -static_surface = mcdc.Surface.TorusZ(A=A, B=B, C=C, R=R, r=r) - -# Test object: moving surface -moving_surface = mcdc.Surface.TorusZ(A=A, B=B, C=C, R=R, r=r) -moving_surface.move(velocities, durations) - -# Create the dummy simulation structure and data -structure_container, data = preparation() -structure = structure_container[0] - -# Get the "compiled" test objects -static_surface = structure["surfaces"][0] -moving_surface = structure["surfaces"][1] - -# Particle object for testing -import mcdc.numba_types as type_ - -particle_container = np.zeros(1, type_.particle_data) -particle = particle_container[0] +static_surface = None +moving_surface = None +data = None +particle_container = None +particle = None # Miscellanies TINY = COINCIDENCE_TOLERANCE * 0.1 # Tiny value within coincidence tolerance @@ -54,6 +40,22 @@ torus_z, ) + +@pytest.fixture(autouse=True) +def setup_geometry_case(compile_surfaces): + global static_surface, moving_surface, data, particle_container, particle + + static_surface_obj = mcdc.Surface.TorusZ(A=A, B=B, C=C, R=R, r=r) + moving_surface_obj = mcdc.Surface.TorusZ(A=A, B=B, C=C, R=R, r=r) + moving_surface_obj.move(velocities, durations) + + static_surface, moving_surface, data = compile_surfaces( + static_surface_obj, moving_surface_obj + ) + particle_container = np.zeros(1, type_.particle_data) + particle = particle_container[0] + + # ===================================================================================== # Torus-Z core functions # ===================================================================================== diff --git a/test/unit/transport/tally/conftest.py b/test/unit/transport/tally/conftest.py index ca58fd7e9..6d1d08703 100644 --- a/test/unit/transport/tally/conftest.py +++ b/test/unit/transport/tally/conftest.py @@ -1,11 +1,6 @@ -import os - import numpy as np import pytest -# Match distribution tests: run Numba-decorated routines in pure Python for unit tests. -os.environ.setdefault("NUMBA_DISABLE_JIT", "1") - import mcdc import mcdc.numba_types as type_ from mcdc.constant import PARTICLE_NEUTRON diff --git a/test/unit/transport/technique/weight_windows.py b/test/unit/transport/technique/test_weight_windows.py similarity index 100% rename from test/unit/transport/technique/weight_windows.py rename to test/unit/transport/technique/test_weight_windows.py diff --git a/test/unit/transport/util/find_bin.py b/test/unit/transport/util/test_find_bin.py similarity index 100% rename from test/unit/transport/util/find_bin.py rename to test/unit/transport/util/test_find_bin.py diff --git a/test/unit/transport/util/linear_interpolation.py b/test/unit/transport/util/test_linear_interpolation.py similarity index 100% rename from test/unit/transport/util/linear_interpolation.py rename to test/unit/transport/util/test_linear_interpolation.py From 4dbd0ecf22fec455c63dfbac1d7a0f7b4ce04fac Mon Sep 17 00:00:00 2001 From: Massimo Larsen Date: Wed, 20 May 2026 13:29:30 -0700 Subject: [PATCH 06/15] Revert "Initial test update commit. Renamed tests, modified distribution tests to run with numba, modified surface tests to run with pytest." This reverts commit f86e9938c4e31118722dc1c066ef86239c33f266. --- pyproject.toml | 1 - test/unit/conftest.py | 22 ----- ...rgy_deposition.py => energy_deposition.py} | 0 test/unit/transport/distributions/conftest.py | 90 ++++++++----------- .../distributions/test_evaporation.py | 23 ++--- .../distributions/test_kalbach_mann.py | 13 ++- .../distributions/test_level_scattering.py | 7 +- .../distributions/test_maxwellian.py | 23 ++--- .../distributions/test_multi_table.py | 23 ++--- .../distributions/test_nbody_correlated.py | 43 +++------ .../test_tabulated_distribution.py | 17 ++-- .../test_tabulated_energy_angle.py | 14 +-- .../transport/geometry/surface/conftest.py | 25 ------ .../surface/{test_plane_x.py => plane_x.py} | 44 +++++---- .../surface/{test_plane_y.py => plane_y.py} | 44 +++++---- .../surface/{test_plane_z.py => plane_z.py} | 44 +++++---- .../surface/{test_torus_z.py => torus_z.py} | 44 +++++---- test/unit/transport/tally/conftest.py | 5 ++ ...st_weight_windows.py => weight_windows.py} | 0 .../util/{test_find_bin.py => find_bin.py} | 0 ...terpolation.py => linear_interpolation.py} | 0 21 files changed, 175 insertions(+), 307 deletions(-) delete mode 100644 test/unit/conftest.py rename test/unit/object_/tally/{test_energy_deposition.py => energy_deposition.py} (100%) delete mode 100644 test/unit/transport/geometry/surface/conftest.py rename test/unit/transport/geometry/surface/{test_plane_x.py => plane_x.py} (96%) rename test/unit/transport/geometry/surface/{test_plane_y.py => plane_y.py} (96%) rename test/unit/transport/geometry/surface/{test_plane_z.py => plane_z.py} (96%) rename test/unit/transport/geometry/surface/{test_torus_z.py => torus_z.py} (97%) rename test/unit/transport/technique/{test_weight_windows.py => weight_windows.py} (100%) rename test/unit/transport/util/{test_find_bin.py => find_bin.py} (100%) rename test/unit/transport/util/{test_linear_interpolation.py => linear_interpolation.py} (100%) diff --git a/pyproject.toml b/pyproject.toml index 2a7c9850c..6e7c63049 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -95,6 +95,5 @@ force-exclude = ''' mcdc/numba_types\.py | mcdc/mcdc_get/ | mcdc/mcdc_set/ - | mcdc/lib/ ) ''' diff --git a/test/unit/conftest.py b/test/unit/conftest.py deleted file mode 100644 index 6ae1f6130..000000000 --- a/test/unit/conftest.py +++ /dev/null @@ -1,22 +0,0 @@ -import os - - -def pytest_addoption(parser): - parser.addoption( - "--mode", - choices=["python", "numba"], - default="python", - help="MCDC execution mode for unit tests.", - ) - - -def pytest_configure(config): - mode = config.getoption("--mode") - if mode == "python": - os.environ["NUMBA_DISABLE_JIT"] = "1" - else: - os.environ.pop("NUMBA_DISABLE_JIT", None) - - -def pytest_report_header(config): - return f"MCDC mode: {config.getoption('--mode')}" diff --git a/test/unit/object_/tally/test_energy_deposition.py b/test/unit/object_/tally/energy_deposition.py similarity index 100% rename from test/unit/object_/tally/test_energy_deposition.py rename to test/unit/object_/tally/energy_deposition.py diff --git a/test/unit/transport/distributions/conftest.py b/test/unit/transport/distributions/conftest.py index 994845027..96ee1b6c9 100644 --- a/test/unit/transport/distributions/conftest.py +++ b/test/unit/transport/distributions/conftest.py @@ -1,66 +1,54 @@ -import numpy as np +import os import pytest -from numba import njit -import mcdc.transport.distribution as dist - -MOCK_RNG_MAX_VALUES = 8 -MOCK_RNG_STATE_DTYPE = np.dtype( - [ - ("idx", np.int64), - ("n_values", np.int64), - ("values", np.float64, (MOCK_RNG_MAX_VALUES,)), - ] -) - - -@njit -def mock_lcg_known_sequence(rng_state): - # Force lcg to pick mock rng values in sequence - state = rng_state[0] - i = state["idx"] - n = state["n_values"] +# Force pure-Python execution for numba in tests. +os.environ.setdefault("NUMBA_DISABLE_JIT", "1") - if i < n: - value = state["values"][i] - else: - value = np.nan - - state["idx"] = i + 1 - return value +import mcdc.transport.distribution as dist -def mock_rng(*values): - if len(values) == 1 and isinstance(values[0], (list, tuple, np.ndarray)): - values = tuple(values[0]) +class MockRNG: + def __init__(self, values): + self._values = list(values) + self._i = 0 - state = np.zeros(1, dtype=MOCK_RNG_STATE_DTYPE) - state[0]["idx"] = 0 - state[0]["n_values"] = len(values) - if values: - state[0]["values"][: len(values)] = np.asarray(values, dtype=np.float64) - return state + def lcg(self, _state_container): + # MCDC uses dist.rng.lcg(state) as its RNG hook. We override it with this + # deterministic sequence so tests can be analytic and reproducible. + if self._i >= len(self._values): + raise IndexError("MockRNG depleted") + value = self._values[self._i] + self._i += 1 + return value @pytest.fixture -def make_distribution_record(): - # Build one typed record from a dict of field values. - def _make(record_dtype, values_dict): - container = np.zeros(1, record_dtype) - record = container[0] - for key, value in values_dict.items(): - record[key] = value - return record +def rng_state(): + """ + Return the minimal RNG state container expected by MCDC. - return _make + MCDC passes a list of per-thread state dicts; tests only need one. + """ + return [dict(rng_seed=0)] @pytest.fixture -def mock_rng_sequence(): - # Replace mcdc lcg with the known sequence +def rng_sequence(): + """ + Temporarily replace dist.rng.lcg with a deterministic sequence. + + Usage: + rng_sequence([0.1, 0.2, 0.3]) + ... code under test that calls dist.rng.lcg(...) + """ original_lcg = dist.rng.lcg - dist.rng.lcg = mock_lcg_known_sequence - # Yield the mock_rng sequence - yield mock_rng - # Reset the lcg. If the test fails, the mock_rng could still be installed and cause later tests to fail + + def _apply(values): + # Swap in a mock LCG implementation that returns the provided values. + rng = MockRNG(values) + dist.rng.lcg = rng.lcg + return rng + + # Yield a callable to the test, then restore the real RNG after the test ends. + yield _apply dist.rng.lcg = original_lcg diff --git a/test/unit/transport/distributions/test_evaporation.py b/test/unit/transport/distributions/test_evaporation.py index 9a5774133..354d131b6 100644 --- a/test/unit/transport/distributions/test_evaporation.py +++ b/test/unit/transport/distributions/test_evaporation.py @@ -1,31 +1,20 @@ import math -import numpy as np -import mcdc.numba_types as type_ import mcdc.transport.distribution as dist from .test_data import make_test_table_data_constant -def test_evaporation_sample(mock_rng_sequence, make_distribution_record): +def test_evaporation_sample(rng_sequence, rng_state): # MCNP Theory & User Manual §2.4.3.5.4.7 (Law 9: Evaporation Spectrum) - table_dict, data = make_test_table_data_constant(1.0) - table = make_distribution_record(type_.table_data, table_dict) - evaporation = make_distribution_record( - type_.evaporation_distribution, - {"nuclear_temperature_ID": 0, "restriction_energy": 0.0}, - ) - data = np.asarray(data, dtype=np.float64) - - simulation_dtype = np.dtype([("table_data", type_.table_data, (1,))]) - simulation_container = np.zeros(1, dtype=simulation_dtype) - simulation = simulation_container[0] - simulation["table_data"][0] = table + table, data = make_test_table_data_constant(1.0) + mcdc = {"table_data": [table]} + evaporation = {"nuclear_temperature_ID": 0, "restriction_energy": 0.0} xi1, xi2 = 0.1, 0.2 - mock_rng = mock_rng_sequence(xi1, xi2) + rng_sequence([xi1, xi2]) - sampled_E = dist.sample_evaporation(2.0, mock_rng, evaporation, simulation, data) + sampled_E = dist.sample_evaporation(2.0, rng_state, evaporation, mcdc, data) # The constant temperature table gives T(E_in) = 1.0 for this test. # With U = 0, the MCNP notation E_in - U becomes 2.0, so # w = (E_in - U) / T(E_in) = 2.0. diff --git a/test/unit/transport/distributions/test_kalbach_mann.py b/test/unit/transport/distributions/test_kalbach_mann.py index 51f912b6f..2d29048c3 100644 --- a/test/unit/transport/distributions/test_kalbach_mann.py +++ b/test/unit/transport/distributions/test_kalbach_mann.py @@ -1,23 +1,20 @@ import math -import numpy as np -import mcdc.numba_types as type_ import mcdc.transport.distribution as dist from .test_data import make_test_kalbach_mann_data -def test_kalbach_mann_sample(mock_rng_sequence, make_distribution_record): +def test_kalbach_mann_sample(rng_sequence, rng_state): # MCNP Theory & User Manual §2.4.3.5.4.11 (Law 44: Kalbach-87 Correlated Energy-angle Scattering) - kalbach_dict, data = make_test_kalbach_mann_data() - kalbach = make_distribution_record(type_.kalbach_mann_distribution, kalbach_dict) - data = np.asarray(data, dtype=np.float64) + kalbach, data = make_test_kalbach_mann_data() # For E_in = 2.0 on the grid [1, 3], the interpolation fraction is r = 0.5. # xi_1 = 0.3 < r, so the second energy table is selected. xi1, xi2, xi3, xi4 = 0.3, 0.1, 0.7, 0.5 - mock_rng = mock_rng_sequence(xi1, xi2, xi3, xi4) - sampled_E, sampled_mu = dist.sample_kalbach_mann(2.0, mock_rng, kalbach, data) + rng_sequence([xi1, xi2, xi3, xi4]) + + sampled_E, sampled_mu = dist.sample_kalbach_mann(2.0, rng_state, kalbach, data) # Law 44 uses the Law 4 energy construction, so with xi_2 = 0.1 in the first bin # of the selected table: diff --git a/test/unit/transport/distributions/test_level_scattering.py b/test/unit/transport/distributions/test_level_scattering.py index 45ca3bf6a..2ecba06e6 100644 --- a/test/unit/transport/distributions/test_level_scattering.py +++ b/test/unit/transport/distributions/test_level_scattering.py @@ -1,18 +1,15 @@ import math -import mcdc.numba_types as type_ import mcdc.transport.distribution as dist -def test_level_scattering_sample(make_distribution_record): +def test_level_scattering_sample(): # MCNP Theory & User Manual §2.4.3.5.4.3 (Law 3: Inelastic Scattering from Nuclear Levels) # The implementation stores the Law 3 relation in the linear form # E_out = C2 * (E_in - C1). # For this test, E_in = 5, C1 = 1, and C2 = 0.5, so # E_out = 0.5 * (5 - 1) = 2. - level = make_distribution_record( - type_.level_scattering_distribution, {"C1": 1.0, "C2": 0.5} - ) + level = {"C1": 1.0, "C2": 0.5} sampled_E = dist.sample_level_scattering(5.0, level) expected_E = 2.0 assert math.isclose(sampled_E, expected_E, rel_tol=0.0, abs_tol=1e-12) diff --git a/test/unit/transport/distributions/test_maxwellian.py b/test/unit/transport/distributions/test_maxwellian.py index 44c85a354..74fa5bec2 100644 --- a/test/unit/transport/distributions/test_maxwellian.py +++ b/test/unit/transport/distributions/test_maxwellian.py @@ -1,32 +1,21 @@ import math -import numpy as np -import mcdc.numba_types as type_ import mcdc.transport.distribution as dist from mcdc.constant import PI from .test_data import make_test_table_data_constant -def test_maxwellian_sample(mock_rng_sequence, make_distribution_record): +def test_maxwellian_sample(rng_sequence, rng_state): # MCNP Theory & User Manual §2.4.3.5.4.6 (Law 7: Simple Maxwell Fission Spectrum) - table_dict, data = make_test_table_data_constant(1.0) - table = make_distribution_record(type_.table_data, table_dict) - maxwellian = make_distribution_record( - type_.maxwellian_distribution, - {"nuclear_temperature_ID": 0, "restriction_energy": 0.0}, - ) - data = np.asarray(data, dtype=np.float64) - - simulation_dtype = np.dtype([("table_data", type_.table_data, (1,))]) - simulation_container = np.zeros(1, dtype=simulation_dtype) - simulation = simulation_container[0] - simulation["table_data"][0] = table + table, data = make_test_table_data_constant(1.0) + mcdc = {"table_data": [table]} + maxwellian = {"nuclear_temperature_ID": 0, "restriction_energy": 0.0} xi1, xi2, xi3 = 0.9, 0.9, 0.0 - mock_rng = mock_rng_sequence(xi1, xi2, xi3) + rng_sequence([xi1, xi2, xi3]) - sampled_E = dist.sample_maxwellian(2.0, mock_rng, maxwellian, simulation, data) + sampled_E = dist.sample_maxwellian(2.0, rng_state, maxwellian, mcdc, data) # MCNP Eq. (2.73): # E_out = -T(E_in) * [xi_1^2 / (xi_1^2 + xi_2^2) * ln(xi_3) + ln(xi_4)] # MCDC uses the equivalent polar-form reduction diff --git a/test/unit/transport/distributions/test_multi_table.py b/test/unit/transport/distributions/test_multi_table.py index bf64c65c1..d1a750758 100644 --- a/test/unit/transport/distributions/test_multi_table.py +++ b/test/unit/transport/distributions/test_multi_table.py @@ -1,31 +1,24 @@ import math -import numpy as np -import mcdc.numba_types as type_ import mcdc.transport.distribution as dist from .test_data import make_test_multi_table_data -def test_multi_table_distribution_sample(mock_rng_sequence, make_distribution_record): +def test_multi_table_distribution_sample(rng_sequence, rng_state): # MCNP Theory & User Manual §2.4.3.5.4.4 (Law 4: Tabular Distribution) - multi_table_dict, data = make_test_multi_table_data() - multi_table = make_distribution_record( - type_.multi_table_distribution, multi_table_dict - ) - data = np.asarray(data, dtype=np.float64) + multi_table, data = make_test_multi_table_data() # For E_in = 2.0 on the grid [1, 3], Eq. (2.62) gives r = 0.5. - # xi_1 = 0.3 < r, so Eq. (2.64) selects l = i + 1, i.e. the second table. - xi1, xi2 = 0.3, 0.2 - mock_rng = mock_rng_sequence(xi1, xi2) + # xi_2 = 0.3 < r, so Eq. (2.64) selects l = i + 1, i.e. the second table. + rng_sequence([0.3, 0.2]) - sampled_E = dist._sample_multi_table(2.0, mock_rng, multi_table, data, scale=True) + sampled_E = dist._sample_multi_table(2.0, rng_state, multi_table, data, scale=True) - # In the selected table, xi_2 = 0.2 falls in the first continuous bin. - # Eq. (2.65) gives E' = E_l,k + (xi_2 - c_l,k) / p_l,k = 100 + 0.2 / 0.01 = 120. + # In the selected table, xi_1 = 0.2 falls in the first continuous bin. + # Eq. (2.65) gives E' = E_l,k + (xi_1 - c_l,k) / p_l,k = 100 + 0.2 / 0.01 = 120. # The test data use constant p within the bin, so the linear-linear form in # Eq. (2.66) reduces to the same result. - E_prime = 100.0 + (xi2 - 0.0) / 0.01 + E_prime = 100.0 + (0.2 - 0.0) / 0.01 # Eq. (2.67) and Eq. (2.68) give the scaled bounds: # E_1 = 10 + 0.5 * (100 - 10) = 55 # E_K = 30 + 0.5 * (300 - 30) = 165 diff --git a/test/unit/transport/distributions/test_nbody_correlated.py b/test/unit/transport/distributions/test_nbody_correlated.py index dc2536783..a35f63361 100644 --- a/test/unit/transport/distributions/test_nbody_correlated.py +++ b/test/unit/transport/distributions/test_nbody_correlated.py @@ -1,48 +1,25 @@ import math -import numpy as np -import mcdc.numba_types as type_ import mcdc.transport.distribution as dist from mcdc.constant import DISTRIBUTION_N_BODY from .test_data import make_test_tabulated_data -def test_nbody_sample_correlated(mock_rng_sequence, make_distribution_record): +def test_nbody_sample_correlated(rng_sequence, rng_state): # MCNP Theory & User Manual §2.4.3.5.4.13 (Law 66: N-body Phase Space Distribution) - table_dict, data = make_test_tabulated_data([2.0, 4.0, 6.0], [0.0, 0.4, 1.0]) - nbody = make_distribution_record(type_.nbody_distribution, table_dict) - distribution = make_distribution_record( - type_.distribution, {"child_type": DISTRIBUTION_N_BODY, "child_ID": 0} - ) - data = np.asarray(data, dtype=np.float64) - - # Numba compiles all correlated-branch field accesses, so this container needs - # the three correlated distribution arrays even though this test uses N-body only. - simulation_dtype = np.dtype( - [ - ("kalbach_mann_distributions", type_.kalbach_mann_distribution, (1,)), - ( - "tabulated_energy_angle_distributions", - type_.tabulated_energy_angle_distribution, - (1,), - ), - ("nbody_distributions", type_.nbody_distribution, (1,)), - ] - ) - simulation_container = np.zeros(1, dtype=simulation_dtype) - simulation = simulation_container[0] - simulation["nbody_distributions"][0] = nbody + table, data = make_test_tabulated_data([2.0, 4.0, 6.0], [0.0, 0.4, 1.0]) + mcdc = {"nbody_distributions": [table]} + distribution = {"child_type": DISTRIBUTION_N_BODY, "child_ID": 0} # First value samples energy, second value samples isotropic cosine. - xi1, xi2 = 0.2, 0.75 - mock_rng = mock_rng_sequence(xi1, xi2) + rng_sequence([0.2, 0.75]) sampled_E, sampled_mu = dist.sample_correlated_distribution( 2.0, distribution, - mock_rng, - simulation, + rng_state, + mcdc, data, ) @@ -50,11 +27,11 @@ def test_nbody_sample_correlated(mock_rng_sequence, make_distribution_record): # samples the cosine isotropically. This test is therefore checking the current # reduced implementation, not reconstructing the full Law 66 rejection sampler # from Eq. (2.103) through Eq. (2.106). - # For the tabulated-energy part, xi_1 = 0.2 gives linear interpolation in the first + # For the tabulated-energy part, xi = 0.2 gives linear interpolation in the first # bin. For the angular part, MCNP Eq. (2.107) gives mu = 2 * xi_10 - 1 for # isotropic center-of-mass sampling. - expected_E = 2.0 + (xi1 - 0.0) * (4.0 - 2.0) / (0.4 - 0.0) - expected_mu = 2.0 * xi2 - 1.0 + expected_E = 2.0 + (0.2 - 0.0) * (4.0 - 2.0) / (0.4 - 0.0) + expected_mu = 2.0 * 0.75 - 1.0 assert math.isclose(sampled_E, expected_E, rel_tol=0.0, abs_tol=1e-12) assert math.isclose(sampled_mu, expected_mu, rel_tol=0.0, abs_tol=1e-12) diff --git a/test/unit/transport/distributions/test_tabulated_distribution.py b/test/unit/transport/distributions/test_tabulated_distribution.py index 93e1d1166..648e1da4f 100644 --- a/test/unit/transport/distributions/test_tabulated_distribution.py +++ b/test/unit/transport/distributions/test_tabulated_distribution.py @@ -1,24 +1,19 @@ import math -import numpy as np -import mcdc.numba_types as type_ import mcdc.transport.distribution as dist from .test_data import make_test_tabulated_data -def test_tabulated_distribution_sample(mock_rng_sequence, make_distribution_record): +def test_tabulated_distribution_sample(rng_sequence, rng_state): # MCNP Theory & User Manual §2.4.3.5.4.4 (Law 4: Tabular Distribution) - table_dict, data = make_test_tabulated_data([1.0, 3.0, 7.0], [0.0, 0.4, 1.0]) - table = make_distribution_record(type_.tabulated_distribution, table_dict) - data = np.asarray(data, dtype=np.float64) - xi1 = 0.2 - mock_rng = mock_rng_sequence(xi1) + table, data = make_test_tabulated_data([1.0, 3.0, 7.0], [0.0, 0.4, 1.0]) + rng_sequence([0.2]) - sampled_E = dist.sample_tabulated(table, mock_rng, data) + sampled_E = dist.sample_tabulated(table, rng_state, data) # This is the single-table inverse-CDF interpolation used by the tabulated sampler: - # xi_1 = 0.2 lies in the first bin, so linear interpolation between + # xi = 0.2 lies in the first bin, so linear interpolation between # (c_0, E_0) = (0.0, 1.0) and (c_1, E_1) = (0.4, 3.0) gives the expected value. - expected_E = 1.0 + (xi1 - 0.0) * (3.0 - 1.0) / (0.4 - 0.0) + expected_E = 1.0 + (0.2 - 0.0) * (3.0 - 1.0) / (0.4 - 0.0) assert math.isclose(sampled_E, expected_E, rel_tol=0.0, abs_tol=1e-12) diff --git a/test/unit/transport/distributions/test_tabulated_energy_angle.py b/test/unit/transport/distributions/test_tabulated_energy_angle.py index 8ed4a5263..8c7d770c7 100644 --- a/test/unit/transport/distributions/test_tabulated_energy_angle.py +++ b/test/unit/transport/distributions/test_tabulated_energy_angle.py @@ -1,25 +1,19 @@ import math -import numpy as np -import mcdc.numba_types as type_ import mcdc.transport.distribution as dist from .test_data import make_test_tabulated_energy_angle_data -def test_tabulated_energy_angle_sample(mock_rng_sequence, make_distribution_record): +def test_tabulated_energy_angle_sample(rng_sequence, rng_state): # MCNP Theory & User Manual §2.4.3.5.4.12 (Law 61: Correlated Energy-angle Scattering) - table_dict, data = make_test_tabulated_energy_angle_data() - table = make_distribution_record( - type_.tabulated_energy_angle_distribution, table_dict - ) - data = np.asarray(data, dtype=np.float64) + table, data = make_test_tabulated_energy_angle_data() xi1, xi2, xi3 = 0.3, 0.1, 0.25 - mock_rng = mock_rng_sequence(xi1, xi2, xi3) + rng_sequence([xi1, xi2, xi3]) sampled_E, sampled_mu = dist.sample_tabulated_energy_angle( - 2.0, mock_rng, table, data + 2.0, rng_state, table, data ) # Law 61 uses the Law 4 energy construction first. diff --git a/test/unit/transport/geometry/surface/conftest.py b/test/unit/transport/geometry/surface/conftest.py deleted file mode 100644 index b26431eda..000000000 --- a/test/unit/transport/geometry/surface/conftest.py +++ /dev/null @@ -1,25 +0,0 @@ -import numpy as np -import pytest - - -@pytest.fixture(autouse=True) -def reset_simulation(): - from mcdc.object_.simulation import simulation - - simulation.__init__() - yield - simulation.__init__() - - -@pytest.fixture -def compile_surfaces(): - def _compile(static_surface_obj, moving_surface_obj): - from mcdc.main import preparation - - structure_container, data = preparation() - structure = structure_container[0] - static_surface = structure["surfaces"][static_surface_obj.ID] - moving_surface = structure["surfaces"][moving_surface_obj.ID] - return static_surface, moving_surface, data - - return _compile diff --git a/test/unit/transport/geometry/surface/test_plane_x.py b/test/unit/transport/geometry/surface/plane_x.py similarity index 96% rename from test/unit/transport/geometry/surface/test_plane_x.py rename to test/unit/transport/geometry/surface/plane_x.py index f52694bdf..20fc2b365 100644 --- a/test/unit/transport/geometry/surface/test_plane_x.py +++ b/test/unit/transport/geometry/surface/plane_x.py @@ -1,6 +1,5 @@ import mcdc import numpy as np -import pytest #### @@ -8,7 +7,7 @@ COINCIDENCE_TOLERANCE, INF, ) -import mcdc.numba_types as type_ +from mcdc.main import preparation # ====================================================================================== # Setup @@ -20,11 +19,26 @@ velocities = np.zeros((3, 3)) velocities[:, 0] = np.array([-1.0, 2.0, -3.0]) -static_surface = None -moving_surface = None -data = None -particle_container = None -particle = None +# Test object: static surface +static_surface = mcdc.Surface.PlaneX(x=X) + +# Test object: moving surface +moving_surface = mcdc.Surface.PlaneX(x=X) +moving_surface.move(velocities, durations) + +# Create the dummy simulation structure and data +structure_container, data = preparation() +structure = structure_container[0] + +# Get the "compiled" test objects +static_surface = structure["surfaces"][0] +moving_surface = structure["surfaces"][1] + +# Particle object for testing +import mcdc.numba_types as type_ + +particle_container = np.zeros(1, type_.particle_data) +particle = particle_container[0] # Miscellanies TINY = COINCIDENCE_TOLERANCE * 0.8 # Tiny value within coincidence tolerance @@ -35,22 +49,6 @@ plane_x, ) - -@pytest.fixture(autouse=True) -def setup_geometry_case(compile_surfaces): - global static_surface, moving_surface, data, particle_container, particle - - static_surface_obj = mcdc.Surface.PlaneX(x=X) - moving_surface_obj = mcdc.Surface.PlaneX(x=X) - moving_surface_obj.move(velocities, durations) - - static_surface, moving_surface, data = compile_surfaces( - static_surface_obj, moving_surface_obj - ) - particle_container = np.zeros(1, type_.particle_data) - particle = particle_container[0] - - # ===================================================================================== # Plane-X core functions # ===================================================================================== diff --git a/test/unit/transport/geometry/surface/test_plane_y.py b/test/unit/transport/geometry/surface/plane_y.py similarity index 96% rename from test/unit/transport/geometry/surface/test_plane_y.py rename to test/unit/transport/geometry/surface/plane_y.py index 7f85629b9..f9e1467d3 100644 --- a/test/unit/transport/geometry/surface/test_plane_y.py +++ b/test/unit/transport/geometry/surface/plane_y.py @@ -1,6 +1,5 @@ import mcdc import numpy as np -import pytest #### @@ -8,7 +7,7 @@ COINCIDENCE_TOLERANCE, INF, ) -import mcdc.numba_types as type_ +from mcdc.main import preparation # ====================================================================================== # Setup @@ -20,11 +19,26 @@ velocities = np.zeros((3, 3)) velocities[:, 1] = np.array([-1.0, 2.0, -3.0]) -static_surface = None -moving_surface = None -data = None -particle_container = None -particle = None +# Test object: static surface +static_surface = mcdc.Surface.PlaneY(y=Y) + +# Test object: moving surface +moving_surface = mcdc.Surface.PlaneY(y=Y) +moving_surface.move(velocities, durations) + +# Create the dummy simulation structure and data +structure_container, data = preparation() +structure = structure_container[0] + +# Get the "compiled" test objects +static_surface = structure["surfaces"][0] +moving_surface = structure["surfaces"][1] + +# Particle object for testing +import mcdc.numba_types as type_ + +particle_container = np.zeros(1, type_.particle_data) +particle = particle_container[0] # Miscellanies TINY = COINCIDENCE_TOLERANCE * 0.8 # Tiny value within coincidence tolerance @@ -35,22 +49,6 @@ plane_y, ) - -@pytest.fixture(autouse=True) -def setup_geometry_case(compile_surfaces): - global static_surface, moving_surface, data, particle_container, particle - - static_surface_obj = mcdc.Surface.PlaneY(y=Y) - moving_surface_obj = mcdc.Surface.PlaneY(y=Y) - moving_surface_obj.move(velocities, durations) - - static_surface, moving_surface, data = compile_surfaces( - static_surface_obj, moving_surface_obj - ) - particle_container = np.zeros(1, type_.particle_data) - particle = particle_container[0] - - # ===================================================================================== # Plane-Y core functions # ===================================================================================== diff --git a/test/unit/transport/geometry/surface/test_plane_z.py b/test/unit/transport/geometry/surface/plane_z.py similarity index 96% rename from test/unit/transport/geometry/surface/test_plane_z.py rename to test/unit/transport/geometry/surface/plane_z.py index 1d68016b2..c4f591ee6 100644 --- a/test/unit/transport/geometry/surface/test_plane_z.py +++ b/test/unit/transport/geometry/surface/plane_z.py @@ -1,6 +1,5 @@ import mcdc import numpy as np -import pytest #### @@ -8,7 +7,7 @@ COINCIDENCE_TOLERANCE, INF, ) -import mcdc.numba_types as type_ +from mcdc.main import preparation # ====================================================================================== # Setup @@ -20,11 +19,26 @@ velocities = np.zeros((3, 3)) velocities[:, 2] = np.array([-1.0, 2.0, -3.0]) -static_surface = None -moving_surface = None -data = None -particle_container = None -particle = None +# Test object: static surface +static_surface = mcdc.Surface.PlaneZ(z=Z) + +# Test object: moving surface +moving_surface = mcdc.Surface.PlaneZ(z=Z) +moving_surface.move(velocities, durations) + +# Create the dummy simulation structure and data +structure_container, data = preparation() +structure = structure_container[0] + +# Get the "compiled" test objects +static_surface = structure["surfaces"][0] +moving_surface = structure["surfaces"][1] + +# Particle object for testing +import mcdc.numba_types as type_ + +particle_container = np.zeros(1, type_.particle_data) +particle = particle_container[0] # Miscellanies TINY = COINCIDENCE_TOLERANCE * 0.8 # Tiny value within coincidence tolerance @@ -35,22 +49,6 @@ plane_z, ) - -@pytest.fixture(autouse=True) -def setup_geometry_case(compile_surfaces): - global static_surface, moving_surface, data, particle_container, particle - - static_surface_obj = mcdc.Surface.PlaneZ(z=Z) - moving_surface_obj = mcdc.Surface.PlaneZ(z=Z) - moving_surface_obj.move(velocities, durations) - - static_surface, moving_surface, data = compile_surfaces( - static_surface_obj, moving_surface_obj - ) - particle_container = np.zeros(1, type_.particle_data) - particle = particle_container[0] - - # ===================================================================================== # Plane-Z core functions # ===================================================================================== diff --git a/test/unit/transport/geometry/surface/test_torus_z.py b/test/unit/transport/geometry/surface/torus_z.py similarity index 97% rename from test/unit/transport/geometry/surface/test_torus_z.py rename to test/unit/transport/geometry/surface/torus_z.py index cbbeba236..c3f83468d 100644 --- a/test/unit/transport/geometry/surface/test_torus_z.py +++ b/test/unit/transport/geometry/surface/torus_z.py @@ -1,7 +1,6 @@ import mcdc import math import numpy as np -import pytest #### @@ -9,7 +8,7 @@ COINCIDENCE_TOLERANCE, INF, ) -import mcdc.numba_types as type_ +from mcdc.main import preparation # ====================================================================================== # Setup @@ -25,11 +24,26 @@ velocities = np.zeros((3, 3)) velocities[:, 0] = np.array([-1.0, 2.0, -3.0]) -static_surface = None -moving_surface = None -data = None -particle_container = None -particle = None +# Test object: static surface +static_surface = mcdc.Surface.TorusZ(A=A, B=B, C=C, R=R, r=r) + +# Test object: moving surface +moving_surface = mcdc.Surface.TorusZ(A=A, B=B, C=C, R=R, r=r) +moving_surface.move(velocities, durations) + +# Create the dummy simulation structure and data +structure_container, data = preparation() +structure = structure_container[0] + +# Get the "compiled" test objects +static_surface = structure["surfaces"][0] +moving_surface = structure["surfaces"][1] + +# Particle object for testing +import mcdc.numba_types as type_ + +particle_container = np.zeros(1, type_.particle_data) +particle = particle_container[0] # Miscellanies TINY = COINCIDENCE_TOLERANCE * 0.1 # Tiny value within coincidence tolerance @@ -40,22 +54,6 @@ torus_z, ) - -@pytest.fixture(autouse=True) -def setup_geometry_case(compile_surfaces): - global static_surface, moving_surface, data, particle_container, particle - - static_surface_obj = mcdc.Surface.TorusZ(A=A, B=B, C=C, R=R, r=r) - moving_surface_obj = mcdc.Surface.TorusZ(A=A, B=B, C=C, R=R, r=r) - moving_surface_obj.move(velocities, durations) - - static_surface, moving_surface, data = compile_surfaces( - static_surface_obj, moving_surface_obj - ) - particle_container = np.zeros(1, type_.particle_data) - particle = particle_container[0] - - # ===================================================================================== # Torus-Z core functions # ===================================================================================== diff --git a/test/unit/transport/tally/conftest.py b/test/unit/transport/tally/conftest.py index 6d1d08703..ca58fd7e9 100644 --- a/test/unit/transport/tally/conftest.py +++ b/test/unit/transport/tally/conftest.py @@ -1,6 +1,11 @@ +import os + import numpy as np import pytest +# Match distribution tests: run Numba-decorated routines in pure Python for unit tests. +os.environ.setdefault("NUMBA_DISABLE_JIT", "1") + import mcdc import mcdc.numba_types as type_ from mcdc.constant import PARTICLE_NEUTRON diff --git a/test/unit/transport/technique/test_weight_windows.py b/test/unit/transport/technique/weight_windows.py similarity index 100% rename from test/unit/transport/technique/test_weight_windows.py rename to test/unit/transport/technique/weight_windows.py diff --git a/test/unit/transport/util/test_find_bin.py b/test/unit/transport/util/find_bin.py similarity index 100% rename from test/unit/transport/util/test_find_bin.py rename to test/unit/transport/util/find_bin.py diff --git a/test/unit/transport/util/test_linear_interpolation.py b/test/unit/transport/util/linear_interpolation.py similarity index 100% rename from test/unit/transport/util/test_linear_interpolation.py rename to test/unit/transport/util/linear_interpolation.py From 1651c0e175b9a700c9e8b709314ac569240792bc Mon Sep 17 00:00:00 2001 From: Massimo Larsen Date: Thu, 21 May 2026 12:47:28 -0700 Subject: [PATCH 07/15] Fixed formatting with newer black --- mcdc/code_factory/numba_objects_generator.py | 6 ++---- 1 file changed, 2 insertions(+), 4 deletions(-) diff --git a/mcdc/code_factory/numba_objects_generator.py b/mcdc/code_factory/numba_objects_generator.py index 1ee8d7306..c466b92d6 100644 --- a/mcdc/code_factory/numba_objects_generator.py +++ b/mcdc/code_factory/numba_objects_generator.py @@ -831,10 +831,8 @@ def align(field_list): pad_id = 0 for field in field_list: if len(field) > 3: - print_error( - "Unexpected struct field specification. Specifications \ - usually only consist of 3 or fewer members" - ) + print_error("Unexpected struct field specification. Specifications \ + usually only consist of 3 or fewer members") multiplier = 1 if len(field) == 3: field = (field[0], field[1], fixup_dims(field[2])) From f58e2944984e76f3cedeff359761822c05acc608 Mon Sep 17 00:00:00 2001 From: Massimo Larsen Date: Fri, 22 May 2026 09:27:18 -0700 Subject: [PATCH 08/15] Updated docs --- docs/source/examples/sphere_in_cube.rst | 2 ++ docs/source/user/first_mcdc.rst | 13 ++++++++++++- 2 files changed, 14 insertions(+), 1 deletion(-) diff --git a/docs/source/examples/sphere_in_cube.rst b/docs/source/examples/sphere_in_cube.rst index b261e8d10..1dca22ee1 100644 --- a/docs/source/examples/sphere_in_cube.rst +++ b/docs/source/examples/sphere_in_cube.rst @@ -126,6 +126,8 @@ Implicit capture is enabled to keep particles alive longer. - Replace the cell filter with a mesh filter to visualise the 3-D flux. - Change the sphere radius or :math:`\nu` to see how fission rate changes. - Add a time grid to the cell tally for time-resolved data. +- Change the cell tally scores to ``["net-current", "current-in", "current-out"]`` + to score current across the sphere cell boundary. Full Input ========== diff --git a/docs/source/user/first_mcdc.rst b/docs/source/user/first_mcdc.rst index 541c224ec..c5067f745 100644 --- a/docs/source/user/first_mcdc.rst +++ b/docs/source/user/first_mcdc.rst @@ -105,7 +105,12 @@ and collision rate. A mesh is created first, then a mesh-filtered ``Tally`` is c Direction bins can also be specified on the tally. Regardless of problem specifics, particles are simulated through all space, direction, and time; the tally definitions are used to indicate in which dimensions a record of particle behavior should be kept. -Available scores include ``"flux"``, ``"density"``, ``"collision"``, ``"capture"``, ``"fission"``, and ``"net-current"``. +Available tracklength scores include ``"flux"``, ``"density"``, ``"collision"``, ``"capture"``, and ``"fission"``. +Current tallies can be attached to either a surface or a cell. +Surface current tallies support ``"net-current"``, while cell current tallies support ``"net-current"``, +``"current-in"``, and ``"current-out"``. +The ``"current-in"`` and ``"current-out"`` scores are positive partial currents; ``"net-current"`` keeps the sign +of the crossing direction. .. code-block:: python3 @@ -117,6 +122,12 @@ Available scores include ``"flux"``, ``"density"``, ``"collision"``, ``"capture" mu=np.linspace(-1.0, 1.0, 32 + 1), ) + # Tally: current crossing a cell boundary + mcdc.Tally( + cell=my_cell, + scores=["net-current", "current-in", "current-out"], + ) + Next we set simulation settings. The only required setting is the number of particles. Settings are configured by assigning attributes on the ``mcdc.settings`` singleton. Additional settings include, for example, the cycles to use for a k-eigenvalue problem From 4164c131ee2fff0184083179fb9449b8dc5abbaa Mon Sep 17 00:00:00 2001 From: Massimo Larsen Date: Mon, 1 Jun 2026 11:47:40 -0700 Subject: [PATCH 09/15] Modified previous TallyCell to be a filter within TallySurface. Updated tests. --- docs/source/examples/sphere_in_cube.rst | 2 +- docs/source/user/first_mcdc.rst | 5 +- mcdc/constant.py | 1 - mcdc/mcdc_get/__init__.py | 2 - mcdc/mcdc_get/cell_tally.py | 3 - mcdc/mcdc_set/__init__.py | 2 - mcdc/mcdc_set/cell_tally.py | 3 - mcdc/numba_types.py | 10 +- mcdc/object_/tally.py | 129 +++++++----------- mcdc/transport/geometry/interface.py | 54 ++------ mcdc/transport/tally/score.py | 82 ++++------- test/unit/transport/tally/test_cell_tally.py | 82 ++++++----- .../tally/test_tally_factory_routing.py | 7 +- 13 files changed, 152 insertions(+), 230 deletions(-) delete mode 100644 mcdc/mcdc_get/cell_tally.py delete mode 100644 mcdc/mcdc_set/cell_tally.py diff --git a/docs/source/examples/sphere_in_cube.rst b/docs/source/examples/sphere_in_cube.rst index 1dca22ee1..2c1033298 100644 --- a/docs/source/examples/sphere_in_cube.rst +++ b/docs/source/examples/sphere_in_cube.rst @@ -126,7 +126,7 @@ Implicit capture is enabled to keep particles alive longer. - Replace the cell filter with a mesh filter to visualise the 3-D flux. - Change the sphere radius or :math:`\nu` to see how fission rate changes. - Add a time grid to the cell tally for time-resolved data. -- Change the cell tally scores to ``["net-current", "current-in", "current-out"]`` +- Change the cell-filtered current scores to ``["net-current", "current-in", "current-out"]`` to score current across the sphere cell boundary. Full Input diff --git a/docs/source/user/first_mcdc.rst b/docs/source/user/first_mcdc.rst index c5067f745..768a0d316 100644 --- a/docs/source/user/first_mcdc.rst +++ b/docs/source/user/first_mcdc.rst @@ -107,8 +107,9 @@ Regardless of problem specifics, particles are simulated through all space, dire the tally definitions are used to indicate in which dimensions a record of particle behavior should be kept. Available tracklength scores include ``"flux"``, ``"density"``, ``"collision"``, ``"capture"``, and ``"fission"``. Current tallies can be attached to either a surface or a cell. -Surface current tallies support ``"net-current"``, while cell current tallies support ``"net-current"``, -``"current-in"``, and ``"current-out"``. +Current scoring uses the surface-crossing estimator in both cases. +For explicit surface filters, supported score is ``"net-current"``. +For cell filters, supported scores are ``"net-current"``, ``"current-in"``, and ``"current-out"``. The ``"current-in"`` and ``"current-out"`` scores are positive partial currents; ``"net-current"`` keeps the sign of the crossing direction. diff --git a/mcdc/constant.py b/mcdc/constant.py index 9bd70766b..c4c62e757 100644 --- a/mcdc/constant.py +++ b/mcdc/constant.py @@ -7,7 +7,6 @@ TALLY_SURFACE = 0 TALLY_COLLISION = 1 TALLY_TRACKLENGTH = 2 -TALLY_CELL = 3 # Tally spatial filter types SPATIAL_FILTER_NONE = 0 diff --git a/mcdc/mcdc_get/__init__.py b/mcdc/mcdc_get/__init__.py index 03dbf2cde..31b31afac 100644 --- a/mcdc/mcdc_get/__init__.py +++ b/mcdc/mcdc_get/__init__.py @@ -102,8 +102,6 @@ import mcdc.mcdc_get.surface_tally as surface_tally -import mcdc.mcdc_get.cell_tally as cell_tally - import mcdc.mcdc_get.collision_tally as collision_tally import mcdc.mcdc_get.tracklength_tally as tracklength_tally diff --git a/mcdc/mcdc_get/cell_tally.py b/mcdc/mcdc_get/cell_tally.py deleted file mode 100644 index fdbf8e750..000000000 --- a/mcdc/mcdc_get/cell_tally.py +++ /dev/null @@ -1,3 +0,0 @@ -# The following is automatically generated by code_factory.py - -from numba import njit diff --git a/mcdc/mcdc_set/__init__.py b/mcdc/mcdc_set/__init__.py index cef3b51c6..e91e83a32 100644 --- a/mcdc/mcdc_set/__init__.py +++ b/mcdc/mcdc_set/__init__.py @@ -102,8 +102,6 @@ import mcdc.mcdc_set.surface_tally as surface_tally -import mcdc.mcdc_set.cell_tally as cell_tally - import mcdc.mcdc_set.collision_tally as collision_tally import mcdc.mcdc_set.tracklength_tally as tracklength_tally diff --git a/mcdc/mcdc_set/cell_tally.py b/mcdc/mcdc_set/cell_tally.py deleted file mode 100644 index fdbf8e750..000000000 --- a/mcdc/mcdc_set/cell_tally.py +++ /dev/null @@ -1,3 +0,0 @@ -# The following is automatically generated by code_factory.py - -from numba import njit diff --git a/mcdc/numba_types.py b/mcdc/numba_types.py index b47342180..5e373cf03 100644 --- a/mcdc/numba_types.py +++ b/mcdc/numba_types.py @@ -674,6 +674,8 @@ ]) surface_tally = into_dtype([ + ('spatial_filter_type', int64), + ('spatial_filter_ID', int64), ('surface_ID', int64), ('filter_surface_bounds', bool), ('has_x_bounds', bool), @@ -689,12 +691,6 @@ ('parent_ID', int64), ]) -cell_tally = into_dtype([ - ('cell_ID', int64), - ('ID', int64), - ('parent_ID', int64), -]) - collision_tally = into_dtype([ ('spatial_filter_type', int64), ('spatial_filter_ID', int64), @@ -844,8 +840,6 @@ def set_simulation(N: dict): ('N_tally', int64), ('surface_tallies', surface_tally, (N['surface_tally'])), ('N_surface_tally', int64), - ('cell_tallies', cell_tally, (N['cell_tally'])), - ('N_cell_tally', int64), ('collision_tallies', collision_tally, (N['collision_tally'])), ('N_collision_tally', int64), ('tracklength_tallies', tracklength_tally, (N['tracklength_tally'])), diff --git a/mcdc/object_/tally.py b/mcdc/object_/tally.py index f4602509c..5c9eeae06 100644 --- a/mcdc/object_/tally.py +++ b/mcdc/object_/tally.py @@ -35,12 +35,12 @@ SCORE_CURRENT_IN, SCORE_CURRENT_OUT, SPATIAL_FILTER_CELL, + SPATIAL_FILTER_SURFACE, SPATIAL_FILTER_MESH, SPATIAL_FILTER_NONE, SURFACE_PLANE_X, SURFACE_PLANE_Y, SURFACE_PLANE_Z, - TALLY_CELL, TALLY_SURFACE, TALLY_COLLISION, TALLY_TRACKLENGTH, @@ -101,7 +101,7 @@ def __new__( x: Iterable[float] | NoneType = None, y: Iterable[float] | NoneType = None, z: Iterable[float] | NoneType = None, - ) -> TallySurface | TallyCell | TallyTracklength | TallyCollision: + ) -> TallySurface | TallyTracklength | TallyCollision: # Determine type and create the tally self based on the provided # spatial filters and scores @@ -124,14 +124,13 @@ def __new__( if surface is None and cell is None: print_error("Current scores need either a surface or a cell tally.") - if surface is not None: - if not set(scores) <= SURFACE_SCORES: - print_error( - "Surface tally currently supports only " - "scores=['net-current']." - ) - return super().__new__(TallySurface) - return super().__new__(TallyCell) + if surface is not None and not set(scores) <= SURFACE_SCORES: + print_error( + "Surface tally currently supports only " "scores=['net-current']." + ) + + # Cell-filtered current tallies share the surface-crossing estimator. + return super().__new__(TallySurface) # Unsupported score for explicit surface selector if surface is not None: @@ -306,8 +305,6 @@ def decode_type(type_): return "Tracklength tally" elif type_ == TALLY_SURFACE: return "Surface tally" - elif type_ == TALLY_CELL: - return "Cell tally" elif type_ == TALLY_COLLISION: return "Collision tally" @@ -333,62 +330,6 @@ def decode_score_type(type_, lower_case=False): return "Current out" if not lower_case else "current-out" -# ====================================================================================== -# Cell tally -# ====================================================================================== - - -class TallyCell(Tally): - # Annotations for Numba mode - label: str = "cell_tally" - # - cell: Cell - - def __init__( - self, - cell: Cell, - name: str = "", - scores: list[str] = ["net-current"], - mu: Iterable[float] | NoneType = None, - azi: Iterable[float] | NoneType = None, - polar_reference: Iterable[float] | NoneType = None, - energy: Iterable[float] | str | NoneType = None, - time: Iterable[float] | NoneType = None, - ): - type_ = TALLY_CELL - super(Tally, self).__init__(type_) - super().__init__( - name, - scores, - mu=mu, - azi=azi, - polar_reference=polar_reference, - energy=energy, - time=time, - ) - - if not set(self.scores) <= { - SCORE_NET_CURRENT, - SCORE_CURRENT_IN, - SCORE_CURRENT_OUT, - }: - print_error( - "Cell tally currently supports only " - "scores=['net-current', 'current-in', 'current-out']." - ) - - # Set cell and attach tally to the cell. - self.cell = cell - cell.tallies.append(self) - - def __repr__(self): - text = super().__repr__() - text += f" - Cell: {self.cell.name}\n" - text += super()._phasespace_filter_text() - text += f" - Bin shape [mu, azi, energy, time, score]: {self.bin_shape} \n" - return text - - # ====================================================================================== # Surface tally # ====================================================================================== @@ -397,7 +338,11 @@ def __repr__(self): class TallySurface(Tally): # Annotations for Numba mode label: str = "surface_tally" + non_numba: list[str] = ["spatial_filter"] # + spatial_filter: Surface | Cell + spatial_filter_type: int + spatial_filter_ID: int surface: Surface filter_surface_bounds: bool has_x_bounds: bool @@ -412,7 +357,8 @@ class TallySurface(Tally): def __init__( self, - surface: Surface, + surface: Surface | NoneType = None, + cell: Cell | NoneType = None, name: str = "", scores: list[str] = ["flux"], mu: Iterable[float] | NoneType = None, @@ -442,9 +388,24 @@ def __init__( "for this tally type." ) - # Set surface and attach tally to the surface - self.surface = surface - surface.tallies.append(self) + if surface is not None and cell is not None: + print_error("Surface tally must specify exactly one of surface or cell.") + if surface is None and cell is None: + print_error("Surface tally requires either a surface or a cell.") + + if surface is not None: + self.spatial_filter = surface + self.spatial_filter_type = SPATIAL_FILTER_SURFACE + self.spatial_filter_ID = surface.ID + self.surface = surface + surface.tallies.append(self) + else: + self.spatial_filter = cell + self.spatial_filter_type = SPATIAL_FILTER_CELL + self.spatial_filter_ID = cell.ID + self.surface = cell.surfaces[0] + for boundary_surface in cell.surfaces: + boundary_surface.tallies.append(self) # Optional bounds for axis-aligned planar surface tallies. self.filter_surface_bounds = False @@ -458,6 +419,13 @@ def __init__( self.z_min = -INF self.z_max = INF + if self.spatial_filter_type == SPATIAL_FILTER_CELL and ( + x is not None or y is not None or z is not None + ): + print_error( + "Cell-filtered surface current tally does not support explicit x/y/z bounds." + ) + if x is not None: x = np.array(x, dtype=float) if len(x) != 2: @@ -493,20 +461,27 @@ def __init__( ) if self.filter_surface_bounds: - if surface.type not in (SURFACE_PLANE_X, SURFACE_PLANE_Y, SURFACE_PLANE_Z): + if self.surface.type not in ( + SURFACE_PLANE_X, + SURFACE_PLANE_Y, + SURFACE_PLANE_Z, + ): print_error( "Bounded surface tally currently supports only PlaneX, PlaneY, and PlaneZ surfaces." ) - if surface.type == SURFACE_PLANE_X and self.has_x_bounds: + if self.surface.type == SURFACE_PLANE_X and self.has_x_bounds: print_error("PlaneX surface tally bounds may only use y and/or z.") - if surface.type == SURFACE_PLANE_Y and self.has_y_bounds: + if self.surface.type == SURFACE_PLANE_Y and self.has_y_bounds: print_error("PlaneY surface tally bounds may only use x and/or z.") - if surface.type == SURFACE_PLANE_Z and self.has_z_bounds: + if self.surface.type == SURFACE_PLANE_Z and self.has_z_bounds: print_error("PlaneZ surface tally bounds may only use x and/or y.") def __repr__(self): text = super().__repr__() - text += f" - Surface: {self.surface.name}\n" + if self.spatial_filter_type == SPATIAL_FILTER_SURFACE: + text += f" - Surface: {self.spatial_filter.name}\n" + elif self.spatial_filter_type == SPATIAL_FILTER_CELL: + text += f" - Cell filter: {self.spatial_filter.name}\n" text += super()._phasespace_filter_text() text += f" - Bin shape [mu, azi, energy, time, score]: {self.bin_shape} \n" return text diff --git a/mcdc/transport/geometry/interface.py b/mcdc/transport/geometry/interface.py index 4b9d049dc..5ddefb509 100644 --- a/mcdc/transport/geometry/interface.py +++ b/mcdc/transport/geometry/interface.py @@ -452,6 +452,9 @@ def surface_crossing(P_arr, simulation, data): elif BC == BC_REFLECTIVE: reflect(P_arr, surface) + pre_cell_ID = -1 + post_cell_ID = -1 + # Score tally for i in range(surface["N_tally"]): tally_ID = int(mcdc_get.surface.tally_IDs(i, surface, data)) @@ -463,22 +466,20 @@ def surface_crossing(P_arr, simulation, data): if not _surface_tally_crossing_in_bounds(P, surface, tally): continue - tally_module.score.surface_tally(P_arr, surface, tally, simulation, data) - - # Score cell net-current tallies tied to this crossed surface. - if BC != BC_REFLECTIVE: - if simulation["N_cell_tally"] > 0: - pre_cell_ID, post_cell_ID = _get_crossing_top_cell_IDs( - P_arr, simulation, data - ) - _score_cell_current_tallies( - P_arr, pre_cell_ID, pre_cell_ID, post_cell_ID, simulation, data - ) - if post_cell_ID != pre_cell_ID: - _score_cell_current_tallies( - P_arr, post_cell_ID, pre_cell_ID, post_cell_ID, simulation, data + # Cell-filtered current tallies need crossing in/out context. + if tally["spatial_filter_type"] == SPATIAL_FILTER_CELL: + if BC == BC_REFLECTIVE: + continue + # only compute this once for each iteration + if pre_cell_ID == -1 and post_cell_ID == -1: + pre_cell_ID, post_cell_ID = _get_crossing_top_cell_IDs( + P_arr, simulation, data ) + tally_module.score.surface_tally( + P_arr, surface, tally, pre_cell_ID, post_cell_ID, simulation, data + ) + # Need to check new cell later? if P["alive"] and not BC == BC_REFLECTIVE: P["cell_ID"] = -1 @@ -514,31 +515,6 @@ def _get_crossing_top_cell_IDs(particle_container, simulation, data): return pre_cell_ID, post_cell_ID -@njit -def _score_cell_current_tallies( - particle_container, - cell_ID, - pre_cell_ID, - post_cell_ID, - simulation, - data, -): - if cell_ID < 0: - return - - cell = simulation["cells"][cell_ID] - for i in range(cell["N_tally"]): - tally_base_ID = int(mcdc_get.cell.tally_IDs(i, cell, data)) - tally_base = simulation["tallies"][tally_base_ID] - if tally_base["child_type"] != TALLY_CELL: - continue - - tally = simulation["cell_tallies"][tally_base["child_ID"]] - tally_module.score.cell_tally( - particle_container, tally, pre_cell_ID, post_cell_ID, simulation, data - ) - - @njit def _surface_tally_crossing_in_bounds(P, surface, tally): surface_type = surface["type"] diff --git a/mcdc/transport/tally/score.py b/mcdc/transport/tally/score.py index ccf6e305e..ffa56b18d 100644 --- a/mcdc/transport/tally/score.py +++ b/mcdc/transport/tally/score.py @@ -27,6 +27,8 @@ SCORE_ENERGY_DEPOSITION, SCORE_CURRENT_IN, SCORE_CURRENT_OUT, + SPATIAL_FILTER_CELL, + SPATIAL_FILTER_SURFACE, SPATIAL_FILTER_MESH, ) from mcdc.transport.geometry.surface import get_normal_component @@ -38,48 +40,9 @@ @njit -def surface_tally(particle_container, surface, tally, simulation, data): - particle = particle_container[0] - tally_base = simulation["tallies"][tally["parent_ID"]] - - # Get filter indices - MG_mode = simulation["settings"]["neutron_multigroup_mode"] - i_mu, i_azi, i_energy, i_time = get_filter_indices( - particle_container, tally_base, data, MG_mode - ) - - # No score if outside non-changing phase-space bins - if i_mu == -1 or i_azi == -1 or i_energy == -1 or i_time == -1: - return - - # Tally index - idx_base = ( - tally_base["bin_offset"] - + i_mu * tally_base["stride_mu"] - + i_azi * tally_base["stride_azi"] - + i_energy * tally_base["stride_energy"] - + i_time * tally_base["stride_time"] - ) - - # Flux - speed = physics.particle_speed(particle_container, simulation, data) - mu = get_normal_component(particle_container, speed, surface, data) - flux = particle["w"] / abs(mu) - - # Score - for i_score in range(tally_base["scores_length"]): - score_type = mcdc_get.tally.scores(i_score, tally_base, data) - score = 0.0 - if score_type == SCORE_NET_CURRENT: - surface = simulation["surfaces"][particle["surface_ID"]] - mu = get_normal_component(particle_container, speed, surface, data) - score = flux * mu - util.atomic_add(data, idx_base + i_score, score) - - -@njit -def cell_tally( +def surface_tally( particle_container, + surface, tally, pre_cell_ID, post_cell_ID, @@ -99,13 +62,6 @@ def cell_tally( if i_mu == -1 or i_azi == -1 or i_energy == -1 or i_time == -1: return - # Incoming/outgoing event flags - target_cell_ID = tally["cell_ID"] - entered = pre_cell_ID != target_cell_ID and post_cell_ID == target_cell_ID - exited = pre_cell_ID == target_cell_ID and post_cell_ID != target_cell_ID - if not entered and not exited: - return - # Tally index idx_base = ( tally_base["bin_offset"] @@ -115,22 +71,42 @@ def cell_tally( + i_time * tally_base["stride_time"] ) + is_cell_filtered = tally["spatial_filter_type"] == SPATIAL_FILTER_CELL + entered = False + exited = False + # check that particle is in correct cell + if is_cell_filtered: + target_cell_ID = tally["spatial_filter_ID"] + entered = pre_cell_ID != target_cell_ID and post_cell_ID == target_cell_ID + exited = pre_cell_ID == target_cell_ID and post_cell_ID != target_cell_ID + if not entered and not exited: + return + + flux = 0.0 + mu = 0.0 + if tally["spatial_filter_type"] == SPATIAL_FILTER_SURFACE: + speed = physics.particle_speed(particle_container, simulation, data) + mu = get_normal_component(particle_container, speed, surface, data) + flux = particle["w"] / abs(mu) + # Score for i_score in range(tally_base["scores_length"]): score_type = mcdc_get.tally.scores(i_score, tally_base, data) score = 0.0 if score_type == SCORE_NET_CURRENT: - if entered: - score = particle["w"] - elif exited: - score = -particle["w"] + if is_cell_filtered: + if entered: + score = particle["w"] + elif exited: + score = -particle["w"] + else: + score = flux * mu elif score_type == SCORE_CURRENT_IN: if entered: score = particle["w"] elif score_type == SCORE_CURRENT_OUT: if exited: score = particle["w"] - if score != 0.0: util.atomic_add(data, idx_base + i_score, score) diff --git a/test/unit/transport/tally/test_cell_tally.py b/test/unit/transport/tally/test_cell_tally.py index e4aeb80d5..56a97452d 100644 --- a/test/unit/transport/tally/test_cell_tally.py +++ b/test/unit/transport/tally/test_cell_tally.py @@ -5,41 +5,45 @@ from mcdc.transport.geometry import interface as geometry_interface -def _bin_value(cell_tally, mcdc_struct, data): - tally_base = mcdc_struct["tallies"][cell_tally["parent_ID"]] +def _bin_value(surface_tally, mcdc_struct, data): + tally_base = mcdc_struct["tallies"][surface_tally["parent_ID"]] return data[tally_base["bin_offset"]] -def _bin_value_score(cell_tally, mcdc_struct, data, score_idx): - tally_base = mcdc_struct["tallies"][cell_tally["parent_ID"]] +def _bin_value_score(surface_tally, mcdc_struct, data, score_idx): + tally_base = mcdc_struct["tallies"][surface_tally["parent_ID"]] return data[tally_base["bin_offset"] + score_idx] -def test_cell_tally_sign_for_incoming_and_outgoing(slab_plane_x, crossing_particle): - cell_tally_obj = mcdc.Tally(cell=slab_plane_x["c_right"], scores=["net-current"]) +def test_surface_cell_filter_current_sign_for_incoming_and_outgoing( + slab_plane_x, crossing_particle +): + current_tally_obj = mcdc.Tally(cell=slab_plane_x["c_right"], scores=["net-current"]) mcdc_container, data = preparation() mcdc_struct = mcdc_container[0] - cell_tally = mcdc_struct["cell_tallies"][cell_tally_obj.child_ID] + current_tally = mcdc_struct["surface_tallies"][current_tally_obj.child_ID] # Left -> right across the shared interior surface: incoming to c_right (+) particle_container = crossing_particle(slab_plane_x["s_mid"].ID, x=0.0, ux=0.5) geometry_interface.surface_crossing(particle_container, mcdc_struct, data) - assert np.isclose(_bin_value(cell_tally, mcdc_struct, data), 2.0) + assert np.isclose(_bin_value(current_tally, mcdc_struct, data), 2.0) # Right -> left across the shared interior surface: outgoing from c_right (-) particle_container = crossing_particle(slab_plane_x["s_mid"].ID, x=0.0, ux=-0.5) geometry_interface.surface_crossing(particle_container, mcdc_struct, data) - assert np.isclose(_bin_value(cell_tally, mcdc_struct, data), 0.0) + assert np.isclose(_bin_value(current_tally, mcdc_struct, data), 0.0) -def test_cell_tally_records_in_and_out_separately(slab_plane_x, crossing_particle): - cell_tally_obj = mcdc.Tally( +def test_surface_cell_filter_current_records_in_and_out_separately( + slab_plane_x, crossing_particle +): + current_tally_obj = mcdc.Tally( cell=slab_plane_x["c_right"], scores=["net-current", "current-in", "current-out"], ) mcdc_container, data = preparation() mcdc_struct = mcdc_container[0] - cell_tally = mcdc_struct["cell_tallies"][cell_tally_obj.child_ID] + current_tally = mcdc_struct["surface_tallies"][current_tally_obj.child_ID] # One incoming and one outgoing crossing. particle_container = crossing_particle(slab_plane_x["s_mid"].ID, x=0.0, ux=0.5) @@ -49,78 +53,84 @@ def test_cell_tally_records_in_and_out_separately(slab_plane_x, crossing_particl # Scores are in requested order: [net-current, current-in, current-out]. # Partial currents are positive; net-current carries the sign. - assert np.isclose(_bin_value_score(cell_tally, mcdc_struct, data, 0), 0.0) - assert np.isclose(_bin_value_score(cell_tally, mcdc_struct, data, 1), 2.0) - assert np.isclose(_bin_value_score(cell_tally, mcdc_struct, data, 2), 2.0) + assert np.isclose(_bin_value_score(current_tally, mcdc_struct, data, 0), 0.0) + assert np.isclose(_bin_value_score(current_tally, mcdc_struct, data, 1), 2.0) + assert np.isclose(_bin_value_score(current_tally, mcdc_struct, data, 2), 2.0) -def test_cell_tally_counts_outgoing_to_vacuum(slab_plane_x, crossing_particle): - cell_tally_obj = mcdc.Tally( +def test_surface_cell_filter_current_counts_outgoing_to_vacuum( + slab_plane_x, crossing_particle +): + current_tally_obj = mcdc.Tally( cell=slab_plane_x["c_right"], scores=["net-current", "current-in", "current-out"], ) mcdc_container, data = preparation() mcdc_struct = mcdc_container[0] - cell_tally = mcdc_struct["cell_tallies"][cell_tally_obj.child_ID] + current_tally = mcdc_struct["surface_tallies"][current_tally_obj.child_ID] # c_right -> vacuum across the right boundary: outgoing (-) particle_container = crossing_particle(slab_plane_x["s_right"].ID, x=1.0, ux=0.5) particle = particle_container[0] geometry_interface.surface_crossing(particle_container, mcdc_struct, data) assert not particle["alive"] - assert np.isclose(_bin_value_score(cell_tally, mcdc_struct, data, 0), -2.0) - assert np.isclose(_bin_value_score(cell_tally, mcdc_struct, data, 1), 0.0) - assert np.isclose(_bin_value_score(cell_tally, mcdc_struct, data, 2), 2.0) + assert np.isclose(_bin_value_score(current_tally, mcdc_struct, data, 0), -2.0) + assert np.isclose(_bin_value_score(current_tally, mcdc_struct, data, 1), 0.0) + assert np.isclose(_bin_value_score(current_tally, mcdc_struct, data, 2), 2.0) -def test_cell_tally_ignores_reflective_crossing(material_mg, crossing_particle): +def test_surface_cell_filter_current_ignores_reflective_crossing( + material_mg, crossing_particle +): s_left = mcdc.Surface.PlaneX(x=-1.0, boundary_condition="vacuum") s_right = mcdc.Surface.PlaneX(x=1.0, boundary_condition="reflective") c_mid = mcdc.Cell(region=+s_left & -s_right, fill=material_mg) - cell_tally_obj = mcdc.Tally( + current_tally_obj = mcdc.Tally( cell=c_mid, scores=["net-current", "current-in", "current-out"], ) mcdc_container, data = preparation() mcdc_struct = mcdc_container[0] - cell_tally = mcdc_struct["cell_tallies"][cell_tally_obj.child_ID] + current_tally = mcdc_struct["surface_tallies"][current_tally_obj.child_ID] particle_container = crossing_particle(s_right.ID, x=1.0, ux=0.5) particle = particle_container[0] geometry_interface.surface_crossing(particle_container, mcdc_struct, data) assert particle["alive"] assert np.isclose(particle["ux"], -0.5) - assert np.isclose(_bin_value_score(cell_tally, mcdc_struct, data, 0), 0.0) - assert np.isclose(_bin_value_score(cell_tally, mcdc_struct, data, 1), 0.0) - assert np.isclose(_bin_value_score(cell_tally, mcdc_struct, data, 2), 0.0) + assert np.isclose(_bin_value_score(current_tally, mcdc_struct, data, 0), 0.0) + assert np.isclose(_bin_value_score(current_tally, mcdc_struct, data, 1), 0.0) + assert np.isclose(_bin_value_score(current_tally, mcdc_struct, data, 2), 0.0) -def test_cell_tally_scores_curved_boundary(material_mg, crossing_particle): +def test_surface_cell_filter_current_scores_curved_boundary( + material_mg, crossing_particle +): s_cyl = mcdc.Surface.CylinderZ(center=(0.0, 0.0), radius=1.0) c_inner = mcdc.Cell(region=-s_cyl, fill=material_mg) mcdc.Cell(region=+s_cyl, fill=material_mg) - cell_tally_obj = mcdc.Tally( + current_tally_obj = mcdc.Tally( cell=c_inner, scores=["net-current", "current-in", "current-out"], ) mcdc_container, data = preparation() mcdc_struct = mcdc_container[0] - cell_tally = mcdc_struct["cell_tallies"][cell_tally_obj.child_ID] + current_tally = mcdc_struct["surface_tallies"][current_tally_obj.child_ID] # Inner -> outer across the cylindrical surface: outgoing from c_inner (-) particle_container = crossing_particle(s_cyl.ID, x=1.0, ux=0.5) geometry_interface.surface_crossing(particle_container, mcdc_struct, data) - assert np.isclose(_bin_value_score(cell_tally, mcdc_struct, data, 0), -2.0) - assert np.isclose(_bin_value_score(cell_tally, mcdc_struct, data, 1), 0.0) - assert np.isclose(_bin_value_score(cell_tally, mcdc_struct, data, 2), 2.0) + assert np.isclose(_bin_value_score(current_tally, mcdc_struct, data, 0), -2.0) + assert np.isclose(_bin_value_score(current_tally, mcdc_struct, data, 1), 0.0) + assert np.isclose(_bin_value_score(current_tally, mcdc_struct, data, 2), 2.0) # Outer -> inner across the same curved surface: incoming to c_inner (+) particle_container = crossing_particle(s_cyl.ID, x=1.0, ux=-0.5) geometry_interface.surface_crossing(particle_container, mcdc_struct, data) - assert np.isclose(_bin_value_score(cell_tally, mcdc_struct, data, 0), 0.0) - assert np.isclose(_bin_value_score(cell_tally, mcdc_struct, data, 1), 2.0) - assert np.isclose(_bin_value_score(cell_tally, mcdc_struct, data, 2), 2.0) + assert np.isclose(_bin_value_score(current_tally, mcdc_struct, data, 0), 0.0) + assert np.isclose(_bin_value_score(current_tally, mcdc_struct, data, 1), 2.0) + assert np.isclose(_bin_value_score(current_tally, mcdc_struct, data, 2), 2.0) diff --git a/test/unit/transport/tally/test_tally_factory_routing.py b/test/unit/transport/tally/test_tally_factory_routing.py index d922d324d..614965fc5 100644 --- a/test/unit/transport/tally/test_tally_factory_routing.py +++ b/test/unit/transport/tally/test_tally_factory_routing.py @@ -3,7 +3,6 @@ import mcdc from mcdc.constant import SPATIAL_FILTER_NONE from mcdc.object_.tally import ( - TallyCell, TallyCollision, TallySurface, TallyTracklength, @@ -12,7 +11,9 @@ def test_tally_factory_routing_surface_vs_tracklength_vs_collision(slab_plane_x): surface_tally = mcdc.Tally(surface=slab_plane_x["s_mid"], scores=["net-current"]) - cell_tally = mcdc.Tally(cell=slab_plane_x["c_right"], scores=["net-current"]) + surface_cell_filter_tally = mcdc.Tally( + cell=slab_plane_x["c_right"], scores=["net-current"] + ) tracklength_tally = mcdc.Tally(scores=["flux"]) mesh = mcdc.MeshUniform( @@ -24,7 +25,7 @@ def test_tally_factory_routing_surface_vs_tracklength_vs_collision(slab_plane_x) collision_tally = mcdc.Tally(mesh=mesh, scores=["energy_deposition"]) assert isinstance(surface_tally, TallySurface) - assert isinstance(cell_tally, TallyCell) + assert isinstance(surface_cell_filter_tally, TallySurface) assert isinstance(tracklength_tally, TallyTracklength) assert isinstance(collision_tally, TallyCollision) assert tracklength_tally.spatial_filter_type == SPATIAL_FILTER_NONE From 7878d8b4a4dcabe235d94982c34ee748f7774418 Mon Sep 17 00:00:00 2001 From: Ilham Variansyah Date: Sun, 7 Jun 2026 10:50:52 +0700 Subject: [PATCH 10/15] update tally score type constants --- mcdc/constant.py | 25 ++++++++++++++----------- 1 file changed, 14 insertions(+), 11 deletions(-) diff --git a/mcdc/constant.py b/mcdc/constant.py index c4c62e757..ee575dc77 100644 --- a/mcdc/constant.py +++ b/mcdc/constant.py @@ -19,22 +19,25 @@ MESH_STRUCTURED = 1 # Tally scores +## Tracklength SCORE_FLUX = 0 SCORE_DENSITY = 1 SCORE_COLLISION = 2 SCORE_CAPTURE = 3 SCORE_FISSION = 4 -SCORE_NET_CURRENT = 5 -SCORE_MU_SQ = 6 -SCORE_TIME_MOMENT_FLUX = 7 -SCORE_SPACE_MOMENT_FLUX = 8 -SCORE_TIME_MOMENT_CURRENT = 9 -SCORE_SPACE_MOMENT_CURRENT = 10 -SCORE_TIME_MOMENT_MU_SQ = 11 -SCORE_SPACE_MOMENT_MU_SQ = 12 -SCORE_ENERGY_DEPOSITION = 13 -SCORE_CURRENT_IN = 14 -SCORE_CURRENT_OUT = 15 +SCORE_MU_SQ = 5 +SCORE_TIME_MOMENT_FLUX = 6 +SCORE_SPACE_MOMENT_FLUX = 7 +SCORE_TIME_MOMENT_CURRENT = 8 +SCORE_SPACE_MOMENT_CURRENT = 9 +SCORE_TIME_MOMENT_MU_SQ = 10 +SCORE_SPACE_MOMENT_MU_SQ = 11 +## Surface-crossing +SCORE_NET_CURRENT = 100 +SCORE_CURRENT_IN = 101 +SCORE_CURRENT_OUT = 102 +## Collision +SCORE_ENERGY_DEPOSITION = 200 # Boundary condition BC_NONE = 0 From ee280eefe0852f28fb10a21b1698627680037657 Mon Sep 17 00:00:00 2001 From: Ilham Variansyah Date: Fri, 12 Jun 2026 05:09:41 +0700 Subject: [PATCH 11/15] Replace Iterable with Sequence --- mcdc/object_/cell.py | 6 ++-- mcdc/object_/mesh.py | 8 ++--- mcdc/object_/source.py | 20 ++++++------ mcdc/object_/surface.py | 20 ++++++------ mcdc/object_/tally.py | 70 ++++++++++++++++++++--------------------- 5 files changed, 62 insertions(+), 62 deletions(-) diff --git a/mcdc/object_/cell.py b/mcdc/object_/cell.py index e2885f642..7c6b6ecd2 100644 --- a/mcdc/object_/cell.py +++ b/mcdc/object_/cell.py @@ -13,7 +13,7 @@ from numpy.typing import NDArray from operator import attrgetter from types import NoneType -from typing import Annotated, Iterable +from typing import Annotated, Sequence from sympy.logic.boolalg import Boolean #### @@ -150,8 +150,8 @@ def __init__( region: Region | NoneType = None, fill: MaterialBase | Universe | Lattice | NoneType = None, name: str = "", - translation: Iterable[float] = [0.0, 0.0, 0.0], - rotation: Iterable[float] = [0.0, 0.0, 0.0], + translation: Sequence[float] = [0.0, 0.0, 0.0], + rotation: Sequence[float] = [0.0, 0.0, 0.0], ): super().__init__() diff --git a/mcdc/object_/mesh.py b/mcdc/object_/mesh.py index 4c915f43e..6307b6ef0 100644 --- a/mcdc/object_/mesh.py +++ b/mcdc/object_/mesh.py @@ -1,4 +1,4 @@ -from typing import Iterable +from typing import Sequence import numpy as np from numpy import float64 @@ -171,9 +171,9 @@ class MeshStructured(MeshBase): def __init__( self, name: str = "", - x: Iterable[float] = [-INF, INF], - y: Iterable[float] = [-INF, INF], - z: Iterable[float] = [-INF, INF], + x: Sequence[float] = [-INF, INF], + y: Sequence[float] = [-INF, INF], + z: Sequence[float] = [-INF, INF], ): type_ = MESH_STRUCTURED super().__init__(type_, name) diff --git a/mcdc/object_/source.py b/mcdc/object_/source.py index 9d8c642e6..697ed7dce 100644 --- a/mcdc/object_/source.py +++ b/mcdc/object_/source.py @@ -3,7 +3,7 @@ from numpy import float64, int64 from numpy.typing import NDArray from types import NoneType -from typing import Annotated, Iterable +from typing import Annotated, Sequence #### @@ -113,21 +113,21 @@ class Source(ObjectNonSingleton): def __init__( self, name: str = "", - position: Iterable[float] | NoneType = None, - x: Iterable[float] | NoneType = None, - y: Iterable[float] | NoneType = None, - z: Iterable[float] | NoneType = None, + position: Sequence[float] | NoneType = None, + x: Sequence[float] | NoneType = None, + y: Sequence[float] | NoneType = None, + z: Sequence[float] | NoneType = None, # - direction: Iterable[float] | NoneType = None, - white_direction: Iterable[float] | NoneType = None, + direction: Sequence[float] | NoneType = None, + white_direction: Sequence[float] | NoneType = None, isotropic: bool | NoneType = None, - polar_cosine: Iterable[float] | NoneType = None, - azimuthal: Iterable[float] | NoneType = None, + polar_cosine: Sequence[float] | NoneType = None, + azimuthal: Sequence[float] | NoneType = None, # energy: float | NDArray[float64] | NoneType = None, energy_group: int | NDArray[int64] | NoneType = None, # - time: float | Iterable[float] = 0.0, + time: float | Sequence[float] = 0.0, # particle_type: str = "neutron", # diff --git a/mcdc/object_/surface.py b/mcdc/object_/surface.py index bbdeaa1d1..19970c54c 100644 --- a/mcdc/object_/surface.py +++ b/mcdc/object_/surface.py @@ -1,4 +1,4 @@ -from typing import Annotated, Iterable +from typing import Annotated, Sequence import numpy as np from numpy import float64 @@ -414,7 +414,7 @@ def Plane( def CylinderX( cls, name: str = "", - center: Iterable[float] = [0.0, 0.0], + center: Sequence[float] = [0.0, 0.0], radius: float = 0.0, boundary_condition: str = "none", ): @@ -458,7 +458,7 @@ def CylinderX( def CylinderY( cls, name: str = "", - center: Iterable[float] = [0.0, 0.0], + center: Sequence[float] = [0.0, 0.0], radius: float = 0.0, boundary_condition: str = "none", ): @@ -502,7 +502,7 @@ def CylinderY( def CylinderZ( cls, name: str = "", - center: Iterable[float] = [0.0, 0.0], + center: Sequence[float] = [0.0, 0.0], radius: float = 0.0, boundary_condition: str = "none", ): @@ -547,8 +547,8 @@ def Cylinder( cls, name: str = "", radius: float = 0.0, - axis: Iterable[float] = [0.0, 0.0, 1.0], - point: Iterable[float] = [0.0, 0.0, 0.0], + axis: Sequence[float] = [0.0, 0.0, 1.0], + point: Sequence[float] = [0.0, 0.0, 0.0], boundary_condition: str = "none", ): """ @@ -603,7 +603,7 @@ def Cylinder( def Sphere( cls, name: str = "", - center: Iterable[float] = [0.0, 0.0, 0.0], + center: Sequence[float] = [0.0, 0.0, 0.0], radius: float = 0.0, boundary_condition: str = "none", ): @@ -649,7 +649,7 @@ def Sphere( def ConeX( cls, name: str = "", - apex: Iterable[float] = [0.0, 0.0, 0.0], + apex: Sequence[float] = [0.0, 0.0, 0.0], t_sq: float = 1.0, boundary_condition: str = "none", ): @@ -693,7 +693,7 @@ def ConeX( def ConeY( cls, name: str = "", - apex: Iterable[float] = [0.0, 0.0, 0.0], + apex: Sequence[float] = [0.0, 0.0, 0.0], t_sq: float = 1.0, boundary_condition: str = "none", ): @@ -736,7 +736,7 @@ def ConeY( def ConeZ( cls, name: str = "", - apex: Iterable[float] = [0.0, 0.0, 0.0], + apex: Sequence[float] = [0.0, 0.0, 0.0], t_sq: float = 1.0, boundary_condition: str = "none", ): diff --git a/mcdc/object_/tally.py b/mcdc/object_/tally.py index 5c9eeae06..c7c5c21eb 100644 --- a/mcdc/object_/tally.py +++ b/mcdc/object_/tally.py @@ -13,7 +13,7 @@ from functools import reduce from numpy import float64 from numpy.typing import NDArray -from typing import Annotated, Iterable +from typing import Annotated, Sequence from types import NoneType #### @@ -93,14 +93,14 @@ def __new__( surface: Surface | NoneType = None, cell: Cell | NoneType = None, mesh: MeshBase | NoneType = None, - mu: Iterable[float] | NoneType = None, - azi: Iterable[float] | NoneType = None, - polar_reference: Iterable[float] | NoneType = None, - energy: Iterable[float] | str | NoneType = None, - time: Iterable[float] | NoneType = None, - x: Iterable[float] | NoneType = None, - y: Iterable[float] | NoneType = None, - z: Iterable[float] | NoneType = None, + mu: Sequence[float] | NoneType = None, + azi: Sequence[float] | NoneType = None, + polar_reference: Sequence[float] | NoneType = None, + energy: Sequence[float] | str | NoneType = None, + time: Sequence[float] | NoneType = None, + x: Sequence[float] | NoneType = None, + y: Sequence[float] | NoneType = None, + z: Sequence[float] | NoneType = None, ) -> TallySurface | TallyTracklength | TallyCollision: # Determine type and create the tally self based on the provided # spatial filters and scores @@ -162,14 +162,14 @@ def __init__( surface: Surface | NoneType = None, cell: Cell | NoneType = None, mesh: MeshBase | NoneType = None, - mu: Iterable[float] | NoneType = None, - azi: Iterable[float] | NoneType = None, - polar_reference: Iterable[float] | NoneType = None, - energy: Iterable[float] | str | NoneType = None, - time: Iterable[float] | NoneType = None, - x: Iterable[float] | NoneType = None, - y: Iterable[float] | NoneType = None, - z: Iterable[float] | NoneType = None, + mu: Sequence[float] | NoneType = None, + azi: Sequence[float] | NoneType = None, + polar_reference: Sequence[float] | NoneType = None, + energy: Sequence[float] | str | NoneType = None, + time: Sequence[float] | NoneType = None, + x: Sequence[float] | NoneType = None, + y: Sequence[float] | NoneType = None, + z: Sequence[float] | NoneType = None, spatial_shape: tuple[int] | NoneType = None, ): # Set name @@ -361,14 +361,14 @@ def __init__( cell: Cell | NoneType = None, name: str = "", scores: list[str] = ["flux"], - mu: Iterable[float] | NoneType = None, - azi: Iterable[float] | NoneType = None, - polar_reference: Iterable[float] | NoneType = None, - energy: Iterable[float] | str | NoneType = None, - time: Iterable[float] | NoneType = None, - x: Iterable[float] | NoneType = None, - y: Iterable[float] | NoneType = None, - z: Iterable[float] | NoneType = None, + mu: Sequence[float] | NoneType = None, + azi: Sequence[float] | NoneType = None, + polar_reference: Sequence[float] | NoneType = None, + energy: Sequence[float] | str | NoneType = None, + time: Sequence[float] | NoneType = None, + x: Sequence[float] | NoneType = None, + y: Sequence[float] | NoneType = None, + z: Sequence[float] | NoneType = None, ): type_ = TALLY_SURFACE super(Tally, self).__init__(type_) @@ -511,11 +511,11 @@ def __init__( mesh: MeshBase | NoneType = None, name: str = "", scores: list[str] = ["energy_deposition"], - mu: Iterable[float] | NoneType = None, - azi: Iterable[float] | NoneType = None, - polar_reference: Iterable[float] | NoneType = None, - energy: Iterable[float] | str | NoneType = None, - time: Iterable[float] | NoneType = None, + mu: Sequence[float] | NoneType = None, + azi: Sequence[float] | NoneType = None, + polar_reference: Sequence[float] | NoneType = None, + energy: Sequence[float] | str | NoneType = None, + time: Sequence[float] | NoneType = None, ): type_ = TALLY_COLLISION spatial_shape = None @@ -619,11 +619,11 @@ def __init__( mesh: MeshBase | NoneType = None, name: str = "", scores: list[str] = ["flux"], - mu: Iterable[float] | NoneType = None, - azi: Iterable[float] | NoneType = None, - polar_reference: Iterable[float] | NoneType = None, - energy: Iterable[float] | str | NoneType = None, - time: Iterable[float] | NoneType = None, + mu: Sequence[float] | NoneType = None, + azi: Sequence[float] | NoneType = None, + polar_reference: Sequence[float] | NoneType = None, + energy: Sequence[float] | str | NoneType = None, + time: Sequence[float] | NoneType = None, ): type_ = TALLY_TRACKLENGTH spatial_shape = None From 08a8873568b7e5cf4cd82da53b67c6d930c49572 Mon Sep 17 00:00:00 2001 From: Ilham Variansyah Date: Fri, 12 Jun 2026 07:33:49 +0700 Subject: [PATCH 12/15] Replace Iterable with Sequence --- mcdc/object_/surface.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/mcdc/object_/surface.py b/mcdc/object_/surface.py index 8b3a59114..f8adf17a6 100644 --- a/mcdc/object_/surface.py +++ b/mcdc/object_/surface.py @@ -1027,8 +1027,8 @@ def TorusZ( def Torus( cls, name: str = "", - center: Iterable[float] = [0.0, 0.0, 0.0], - axis: Iterable[float] = [0.0, 0.0, 1.0], + center: Sequence[float] = [0.0, 0.0, 0.0], + axis: Sequence[float] = [0.0, 0.0, 1.0], R: float = 0.0, r: float = 0.0, boundary_condition: str = "none", From 6f6a40640b0f6acda8b44fc0964148375968024c Mon Sep 17 00:00:00 2001 From: Ilham Variansyah Date: Fri, 12 Jun 2026 08:55:12 +0700 Subject: [PATCH 13/15] Rename TallySurface to TallySurfaceCrossing --- mcdc/object_/surface.py | 4 ++-- mcdc/object_/tally.py | 6 +++--- 2 files changed, 5 insertions(+), 5 deletions(-) diff --git a/mcdc/object_/surface.py b/mcdc/object_/surface.py index f8adf17a6..ececf78be 100644 --- a/mcdc/object_/surface.py +++ b/mcdc/object_/surface.py @@ -31,7 +31,7 @@ ) from mcdc.object_.base import ObjectNonSingleton from mcdc.object_.cell import Region -from mcdc.object_.tally import TallySurface +from mcdc.object_.tally import TallySurfaceCrossing from mcdc.object_.util import move_object from mcdc.print_ import print_error @@ -132,7 +132,7 @@ class Surface(ObjectNonSingleton): move_durations: Annotated[NDArray[float64], ("N_move",)] move_time_grid: Annotated[NDArray[float64], ("N_move_grid",)] move_translations: Annotated[NDArray[float64], ("N_move_grid", 3)] - tallies: list[TallySurface] + tallies: list[TallySurfaceCrossing] def __init__(self, type_, name, boundary_condition): super().__init__() diff --git a/mcdc/object_/tally.py b/mcdc/object_/tally.py index c7c5c21eb..503db322a 100644 --- a/mcdc/object_/tally.py +++ b/mcdc/object_/tally.py @@ -101,7 +101,7 @@ def __new__( x: Sequence[float] | NoneType = None, y: Sequence[float] | NoneType = None, z: Sequence[float] | NoneType = None, - ) -> TallySurface | TallyTracklength | TallyCollision: + ) -> TallySurfaceCrossing | TallyTracklength | TallyCollision: # Determine type and create the tally self based on the provided # spatial filters and scores @@ -130,7 +130,7 @@ def __new__( ) # Cell-filtered current tallies share the surface-crossing estimator. - return super().__new__(TallySurface) + return super().__new__(TallySurfaceCrossing) # Unsupported score for explicit surface selector if surface is not None: @@ -335,7 +335,7 @@ def decode_score_type(type_, lower_case=False): # ====================================================================================== -class TallySurface(Tally): +class TallySurfaceCrossing(Tally): # Annotations for Numba mode label: str = "surface_tally" non_numba: list[str] = ["spatial_filter"] From 43e7a2b38f3494c91335e0147c861ed0dde4e8c3 Mon Sep 17 00:00:00 2001 From: Ilham Variansyah Date: Fri, 12 Jun 2026 08:57:37 +0700 Subject: [PATCH 14/15] Replace TALLY_SURFACE with TALLY_SURFACE_CROSSING --- mcdc/constant.py | 2 +- mcdc/object_/tally.py | 8 ++++---- 2 files changed, 5 insertions(+), 5 deletions(-) diff --git a/mcdc/constant.py b/mcdc/constant.py index 452b49029..b67986e6c 100644 --- a/mcdc/constant.py +++ b/mcdc/constant.py @@ -4,7 +4,7 @@ TALLY = 0 # Tally types -TALLY_SURFACE = 0 +TALLY_SURFACE_CROSSING = 0 TALLY_COLLISION = 1 TALLY_TRACKLENGTH = 2 diff --git a/mcdc/object_/tally.py b/mcdc/object_/tally.py index 503db322a..d18e49dba 100644 --- a/mcdc/object_/tally.py +++ b/mcdc/object_/tally.py @@ -41,7 +41,7 @@ SURFACE_PLANE_X, SURFACE_PLANE_Y, SURFACE_PLANE_Z, - TALLY_SURFACE, + TALLY_SURFACE_CROSSING, TALLY_COLLISION, TALLY_TRACKLENGTH, ) @@ -303,8 +303,8 @@ def __repr__(self): def decode_type(type_): if type_ == TALLY_TRACKLENGTH: return "Tracklength tally" - elif type_ == TALLY_SURFACE: - return "Surface tally" + elif type_ == TALLY_SURFACE_CROSSING: + return "Surface crossing tally" elif type_ == TALLY_COLLISION: return "Collision tally" @@ -370,7 +370,7 @@ def __init__( y: Sequence[float] | NoneType = None, z: Sequence[float] | NoneType = None, ): - type_ = TALLY_SURFACE + type_ = TALLY_SURFACE_CROSSING super(Tally, self).__init__(type_) super().__init__( name, From 47ce31c656d6bc50e6bd2ade84766e57025e3119 Mon Sep 17 00:00:00 2001 From: Ilham Variansyah Date: Fri, 12 Jun 2026 08:59:12 +0700 Subject: [PATCH 15/15] Update unit test for surface crossing tally --- test/unit/transport/tally/test_tally_factory_routing.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/test/unit/transport/tally/test_tally_factory_routing.py b/test/unit/transport/tally/test_tally_factory_routing.py index 614965fc5..ee03f1bfc 100644 --- a/test/unit/transport/tally/test_tally_factory_routing.py +++ b/test/unit/transport/tally/test_tally_factory_routing.py @@ -4,7 +4,7 @@ from mcdc.constant import SPATIAL_FILTER_NONE from mcdc.object_.tally import ( TallyCollision, - TallySurface, + TallySurfaceCrossing, TallyTracklength, ) @@ -24,8 +24,8 @@ def test_tally_factory_routing_surface_vs_tracklength_vs_collision(slab_plane_x) ) collision_tally = mcdc.Tally(mesh=mesh, scores=["energy_deposition"]) - assert isinstance(surface_tally, TallySurface) - assert isinstance(surface_cell_filter_tally, TallySurface) + assert isinstance(surface_tally, TallySurfaceCrossing) + assert isinstance(surface_cell_filter_tally, TallySurfaceCrossing) assert isinstance(tracklength_tally, TallyTracklength) assert isinstance(collision_tally, TallyCollision) assert tracklength_tally.spatial_filter_type == SPATIAL_FILTER_NONE