Add current tallies for bounded surfaces and cells

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
 ==========

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

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

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),
 ])

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__()
 

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)

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",
         #

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",