diff --git a/docs/source/examples/sphere_in_cube.rst b/docs/source/examples/sphere_in_cube.rst index b261e8d10..2c1033298 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-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 541c224ec..768a0d316 100644 --- a/docs/source/user/first_mcdc.rst +++ b/docs/source/user/first_mcdc.rst @@ -105,7 +105,13 @@ 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. +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. .. code-block:: python3 @@ -117,6 +123,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 diff --git a/mcdc/constant.py b/mcdc/constant.py index a4a03218c..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 @@ -19,20 +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_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 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/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 1caa25648..fad9e778c 100644 --- a/mcdc/numba_types.py +++ b/mcdc/numba_types.py @@ -677,7 +677,19 @@ ]) surface_tally = into_dtype([ + ('spatial_filter_type', int64), + ('spatial_filter_ID', int64), ('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_/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 f02e3ff57..ececf78be 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 @@ -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__() @@ -452,7 +452,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", ): @@ -498,7 +498,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", ): @@ -544,7 +544,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", ): @@ -592,8 +592,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", ): """ @@ -651,7 +651,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", ): @@ -699,7 +699,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", ): @@ -746,7 +746,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", ): @@ -792,7 +792,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", ): @@ -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", diff --git a/mcdc/object_/tally.py b/mcdc/object_/tally.py index f55632bbd..d18e49dba 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 #### @@ -32,10 +32,16 @@ SCORE_FISSION, SCORE_NET_CURRENT, SCORE_ENERGY_DEPOSITION, + SCORE_CURRENT_IN, + SCORE_CURRENT_OUT, SPATIAL_FILTER_CELL, + SPATIAL_FILTER_SURFACE, SPATIAL_FILTER_MESH, SPATIAL_FILTER_NONE, - TALLY_SURFACE, + SURFACE_PLANE_X, + SURFACE_PLANE_Y, + SURFACE_PLANE_Z, + TALLY_SURFACE_CROSSING, TALLY_COLLISION, TALLY_TRACKLENGTH, ) @@ -44,9 +50,10 @@ from mcdc.object_.simulation import simulation from mcdc.print_ import print_1d_array, print_error -SURFACE_SCORES = set(["net-current"]) -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): @@ -86,24 +93,52 @@ 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, - ) -> TallySurface | TallyTracklength | TallyCollision: + 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, + ) -> TallySurfaceCrossing | 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_current_score = any(score in CELL_CURRENT_SCORES for score in scores) + + # Surface/cell current tally + if has_current_score: for score in scores: - if not score in SURFACE_SCORES: + if score not in CELL_CURRENT_SCORES: print_error( - f"Scoring '{score}' with surface tally is not supported. " - f"Supported surface tally scores: {SURFACE_SCORES}." + "Cannot mix current scores with non-current scores " + f"in one tally. Current scores: {CELL_CURRENT_SCORES}." ) - return super().__new__(TallySurface) + + if surface is not None and cell is not None: + print_error( + "Current tally must specify exactly one of surface or cell." + ) + + 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 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__(TallySurfaceCrossing) + + # 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: @@ -113,11 +148,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." @@ -132,11 +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, + 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 @@ -160,6 +193,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: @@ -266,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" @@ -287,6 +324,10 @@ 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" # ====================================================================================== @@ -294,24 +335,42 @@ 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"] # + spatial_filter: Surface | Cell + spatial_filter_type: int + spatial_filter_ID: int 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, - surface: Surface, + surface: Surface | NoneType = None, + 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, + 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 + type_ = TALLY_SURFACE_CROSSING super(Tally, self).__init__(type_) super().__init__( name, @@ -329,13 +388,100 @@ 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 + 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 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: + 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 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 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 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 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 @@ -365,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 @@ -473,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 diff --git a/mcdc/transport/geometry/interface.py b/mcdc/transport/geometry/interface.py index 0de3f79aa..5ddefb509 100644 --- a/mcdc/transport/geometry/interface.py +++ b/mcdc/transport/geometry/interface.py @@ -441,20 +441,44 @@ 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"]] + surface = simulation["surfaces"][crossed_surface_ID] BC = surface["boundary_condition"] + + # Apply BC if BC == BC_VACUUM: P["alive"] = False 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)) tally = simulation["surface_tallies"][tally_ID] - tally_module.score.surface_tally(P_arr, surface, tally, simulation, data) + + # 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 + + # 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: @@ -462,6 +486,70 @@ def surface_crossing(P_arr, simulation, data): 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 _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/mcdc/transport/tally/score.py b/mcdc/transport/tally/score.py index 924c03386..ffa56b18d 100644 --- a/mcdc/transport/tally/score.py +++ b/mcdc/transport/tally/score.py @@ -25,6 +25,10 @@ SCORE_FISSION, SCORE_NET_CURRENT, 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 @@ -36,7 +40,15 @@ @njit -def surface_tally(particle_container, surface, tally, simulation, data): +def surface_tally( + particle_container, + surface, + tally, + pre_cell_ID, + post_cell_ID, + simulation, + data, +): particle = particle_container[0] tally_base = simulation["tallies"][tally["parent_ID"]] @@ -59,20 +71,44 @@ def surface_tally(particle_container, surface, tally, simulation, data): + 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) + 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: - 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) + 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/__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..ca58fd7e9 --- /dev/null +++ b/test/unit/transport/tally/conftest.py @@ -0,0 +1,123 @@ +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 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. + 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_cell_tally.py b/test/unit/transport/tally/test_cell_tally.py new file mode 100644 index 000000000..56a97452d --- /dev/null +++ b/test/unit/transport/tally/test_cell_tally.py @@ -0,0 +1,136 @@ +import numpy as np + +import mcdc +from mcdc.main import preparation +from mcdc.transport.geometry import interface as geometry_interface + + +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(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_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] + 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(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(current_tally, mcdc_struct, data), 0.0) + + +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] + 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) + geometry_interface.surface_crossing(particle_container, mcdc_struct, data) + 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]. + # Partial currents are positive; net-current carries the sign. + 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_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] + 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(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_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) + current_tally_obj = mcdc.Tally( + cell=c_mid, + scores=["net-current", "current-in", "current-out"], + ) + + mcdc_container, data = preparation() + mcdc_struct = mcdc_container[0] + 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(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_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) + current_tally_obj = mcdc.Tally( + cell=c_inner, + scores=["net-current", "current-in", "current-out"], + ) + + mcdc_container, data = preparation() + mcdc_struct = mcdc_container[0] + 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(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(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_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..03e59c51e --- /dev/null +++ b/test/unit/transport/tally/test_nonbounded_surface_tally.py @@ -0,0 +1,99 @@ +import numpy as np +import pytest + +import mcdc +from mcdc.main import preparation +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) + + +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_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..ee03f1bfc --- /dev/null +++ b/test/unit/transport/tally/test_tally_factory_routing.py @@ -0,0 +1,61 @@ +import pytest + +import mcdc +from mcdc.constant import SPATIAL_FILTER_NONE +from mcdc.object_.tally import ( + TallyCollision, + TallySurfaceCrossing, + 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"]) + surface_cell_filter_tally = mcdc.Tally( + cell=slab_plane_x["c_right"], 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, 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 + + +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 + + +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 "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): + 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 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