Merge pull request #243 from ethan-lame/dev

Cell Tallies

examples/fixed_source/cooper1/input.py

@@ -0,0 +1,63 @@
+import numpy as np
+import mcdc
+
+
+# =============================================================================
+# Set model
+# =============================================================================
+# A problem based on Problem 1 of [Coper NSE 2001]
+# https://ans.tandfonline.com/action/showCitFormats?doi=10.13182/NSE00-34
+
+# Set materials
+SigmaT = 0.01
+c = 0.8
+m_duct = mcdc.material(capture=np.array([SigmaT]), scatter=np.array([[SigmaT * c]]))
+SigmaT = 1.0
+m_room = mcdc.material(capture=np.array([SigmaT]), scatter=np.array([[SigmaT * c]]))
+
+# Set surfaces
+sx1 = mcdc.surface("plane-x", x=0.0, bc="reflective")
+sx2 = mcdc.surface("plane-x", x=2.0)
+sx3 = mcdc.surface("plane-x", x=2.4)
+sx4 = mcdc.surface("plane-x", x=4.0, bc="vacuum")
+sy1 = mcdc.surface("plane-y", y=0.0, bc="reflective")
+sy2 = mcdc.surface("plane-y", y=1.0)
+sy3 = mcdc.surface("plane-y", y=2.6)
+sy4 = mcdc.surface("plane-y", y=3)
+sy5 = mcdc.surface("plane-y", y=4.0, bc="vacuum")
+
+# Set cells
+# Room
+mcdc.cell(+sx1 & -sx2 & +sy1 & -sy5, m_room)
+mcdc.cell(+sx1 & -sx4 & +sy4 & -sy5, m_room)
+mcdc.cell(+sx2 & -sx3 & +sy1 & -sy2, m_room)
+mcdc.cell(+sx3 & -sx4 & +sy1 & -sy3, m_room)
+
+# Duct Channel
+channel_1 = +sx2 & -sx3 & +sy2 & -sy3
+channel_2 = +sx2 & -sx4 & +sy3 & -sy4
+mcdc.cell(channel_1 | channel_2, m_duct)
+
+# =============================================================================
+# Set source
+# =============================================================================
+# Uniform isotropic source throughout the domain
+
+mcdc.source(x=[0.0, 1.0], y=[0.0, 1.0], isotropic=True)
+
+# =============================================================================
+# Set tally, setting, and run mcdc
+# =============================================================================
+
+mcdc.tally.mesh_tally(
+    scores=["flux"],
+    x=np.linspace(0.0, 4.0, 40),
+    y=np.linspace(0.0, 4.0, 40),
+)
+
+# Setting
+mcdc.setting(N_particle=50)
+mcdc.implicit_capture()
+
+# Run
+mcdc.run()

examples/fixed_source/cooper1/process.py

@@ -0,0 +1,38 @@
+import matplotlib.pyplot as plt
+import h5py
+import numpy as np
+
+# Load result
+with h5py.File("output.h5", "r") as f:
+    tally = f["tallies/mesh_tally_0"]
+    x = tally["grid/x"][:]
+    dx = [x[1:] - x[:-1]][-1]
+    x_mid = 0.5 * (x[:-1] + x[1:])
+
+    phi = tally["flux/mean"][:]
+    phi_sd = tally["flux/sdev"][:]
+
+
+# Plot result
+X, Y = np.meshgrid(x_mid, x_mid)
+Z = phi
+Z = np.log10(np.abs(Z))
+fig, ax = plt.subplots(subplot_kw={"projection": "3d"})
+ax.plot_surface(Y, X, Z, edgecolor="k", color="white")
+ax.set_xlabel("x")
+ax.set_ylabel("y")
+ax.zaxis.set_rotate_label(False)
+ax.set_zlabel("Log10 of scalar flux", rotation=90)
+ax.view_init(elev=18, azim=38)
+plt.show()
+
+Z = phi_sd / phi
+Z = np.log10(np.abs(Z))
+fig, ax = plt.subplots(subplot_kw={"projection": "3d"})
+ax.plot_surface(Y, X, Z, edgecolor="k", color="white")
+ax.set_xlabel("x")
+ax.set_ylabel("y")
+ax.zaxis.set_rotate_label(False)
+ax.set_zlabel("Log10 of scalar flux rel. stdev.", rotation=90)
+ax.view_init(elev=18, azim=38)
+plt.show()

examples/fixed_source/cooper2/input.py

@@ -41,13 +41,14 @@
 # Set tally, setting, and run mcdc
 # =============================================================================
 
-# Tally: cell-average and cell-edge angular fluxes and currents
 mcdc.tally.mesh_tally(
-    scores=["flux"], x=np.linspace(0.0, 4.0, 40), y=np.linspace(0.0, 4.0, 40)
+    scores=["flux"],
+    x=np.linspace(0.0, 4.0, 40),
+    y=np.linspace(0.0, 4.0, 40),
 )
 
 # Setting
-mcdc.setting(N_particle=1e2)
+mcdc.setting(N_particle=50)
 mcdc.implicit_capture()
 
 # Run

examples/fixed_source/cooper_combo/input.py

@@ -0,0 +1,68 @@
+import numpy as np
+import mcdc
+
+
+# =============================================================================
+# Set model
+# =============================================================================
+# A problem based on a combination of Problems 1 & 2 of [Coper NSE 2001]
+# https://ans.tandfonline.com/action/showCitFormats?doi=10.13182/NSE00-34
+
+# Set materials
+SigmaT = 5.0
+c = 0.8
+m_barrier = mcdc.material(capture=np.array([SigmaT]), scatter=np.array([[SigmaT * c]]))
+SigmaT = 1.0
+m_room = mcdc.material(capture=np.array([SigmaT]), scatter=np.array([[SigmaT * c]]))
+SigmaT = 0.01
+m_duct = mcdc.material(capture=np.array([SigmaT]), scatter=np.array([[SigmaT * c]]))
+
+
+# Set surfaces
+sx1 = mcdc.surface("plane-x", x=0.0, bc="reflective")
+sx2 = mcdc.surface("plane-x", x=2.0)
+sx3 = mcdc.surface("plane-x", x=2.4)
+sx4 = mcdc.surface("plane-x", x=4.0, bc="vacuum")
+sy1 = mcdc.surface("plane-y", y=0.0, bc="reflective")
+sy2 = mcdc.surface("plane-y", y=1.0)
+sy3 = mcdc.surface("plane-y", y=2.6)
+sy4 = mcdc.surface("plane-y", y=3)
+sy5 = mcdc.surface("plane-y", y=4.0, bc="vacuum")
+
+# Set cells
+# Room
+mcdc.cell(+sx1 & -sx2 & +sy1 & -sy5, m_room)
+mcdc.cell(+sx1 & -sx4 & +sy4 & -sy5, m_room)
+mcdc.cell(+sx3 & -sx4 & +sy1 & -sy3, m_room)
+
+# Barrier
+mcdc.cell(+sx2 & -sx3 & +sy1 & -sy2, m_barrier)
+
+# Duct Channel
+channel_1 = +sx2 & -sx3 & +sy2 & -sy3
+channel_2 = +sx2 & -sx4 & +sy3 & -sy4
+mcdc.cell(channel_1 | channel_2, m_duct)
+
+# =============================================================================
+# Set source
+# =============================================================================
+# Uniform isotropic source throughout the domain
+
+mcdc.source(x=[0.0, 1.0], y=[0.0, 1.0], isotropic=True)
+
+# =============================================================================
+# Set tally, setting, and run mcdc
+# =============================================================================
+
+mcdc.tally.mesh_tally(
+    scores=["flux"],
+    x=np.linspace(0.0, 4.0, 40),
+    y=np.linspace(0.0, 4.0, 40),
+)
+
+# Setting
+mcdc.setting(N_particle=50)
+mcdc.implicit_capture()
+
+# Run
+mcdc.run()

examples/fixed_source/cooper_combo/process.py

@@ -0,0 +1,39 @@
+import matplotlib.pyplot as plt
+import h5py
+import numpy as np
+
+
+# Load result
+with h5py.File("output.h5", "r") as f:
+    tally = f["tallies/mesh_tally_0"]
+    x = tally["grid/x"][:]
+    dx = [x[1:] - x[:-1]][-1]
+    x_mid = 0.5 * (x[:-1] + x[1:])
+
+    phi = tally["flux/mean"][:]
+    phi_sd = tally["flux/sdev"][:]
+
+
+# Plot result
+X, Y = np.meshgrid(x_mid, x_mid)
+Z = phi
+Z = np.log10(np.abs(Z))
+fig, ax = plt.subplots(subplot_kw={"projection": "3d"})
+ax.plot_surface(Y, X, Z, edgecolor="k", color="white")
+ax.set_xlabel("x")
+ax.set_ylabel("y")
+ax.zaxis.set_rotate_label(False)
+ax.set_zlabel("Log10 of scalar flux", rotation=90)
+ax.view_init(elev=18, azim=38)
+plt.show()
+
+Z = phi_sd / phi
+Z = np.log10(np.abs(Z))
+fig, ax = plt.subplots(subplot_kw={"projection": "3d"})
+ax.plot_surface(Y, X, Z, edgecolor="k", color="white")
+ax.set_xlabel("x")
+ax.set_ylabel("y")
+ax.zaxis.set_rotate_label(False)
+ax.set_zlabel("Log10 of scalar flux rel. stdev.", rotation=90)
+ax.view_init(elev=18, azim=38)
+plt.show()

examples/fixed_source/kobayashi3-TD/input.py

@@ -30,17 +30,17 @@
 
 # Set cells
 # Source
-mcdc.cell(+sx1 & -sx2 & +sy1 & -sy2 & +sz1 & -sz2, m)
+source_cell = mcdc.cell(+sx1 & -sx2 & +sy1 & -sy2 & +sz1 & -sz2, m)
 # Voids
 channel_1 = +sx1 & -sx2 & +sy2 & -sy3 & +sz1 & -sz2
 channel_2 = +sx1 & -sx3 & +sy3 & -sy4 & +sz1 & -sz2
 channel_3 = +sx3 & -sx4 & +sy3 & -sy4 & +sz1 & -sz3
 channel_4 = +sx3 & -sx4 & +sy3 & -sy5 & +sz3 & -sz4
 void_channel = channel_1 | channel_2 | channel_3 | channel_4
-mcdc.cell(void_channel, m_void)
+void_cell = mcdc.cell(void_channel, m_void)
 # Shield
 box = +sx1 & -sx5 & +sy1 & -sy5 & +sz1 & -sz5
-mcdc.cell(box & ~void_channel, m)
+shield_cell = mcdc.cell(box & ~void_channel, m)
 
 # =============================================================================
 # Set source
@@ -58,13 +58,17 @@
 # Tally: z-integrated flux (X-Y section view)
 mcdc.tally.mesh_tally(
     scores=["flux"],
-    x=np.linspace(0.0, 60.0, 61),
-    y=np.linspace(0.0, 100.0, 101),
+    x=np.linspace(0.0, 60.0, 31),
+    y=np.linspace(0.0, 100.0, 51),
     t=np.linspace(0.0, 200.0, 21),
 )
 
+mcdc.tally.cell_tally(source_cell, scores=["flux"])
+mcdc.tally.cell_tally(void_cell, scores=["flux"])
+mcdc.tally.cell_tally(shield_cell, scores=["flux"])
+
 # Setting
-mcdc.setting(N_particle=1e2)
+mcdc.setting(N_particle=80)
 
 # Run
 mcdc.run()

examples/fixed_source/kobayashi3-TD/process.py

@@ -25,6 +25,12 @@
     phi = flux["mean"][:]
     phi_sd = flux["sdev"][:]
 
+    for i in range(len(f["input_deck"]["cell_tallies"])):
+        flux_score = f[f"tallies/cell_tally_{i}/flux"]
+        print(
+            f'cell {i+1} mean = {flux_score["mean"][()]}, sdev = {flux_score["sdev"][()]}'
+        )
+
 fig, ax = plt.subplots()
 cax = ax.pcolormesh(X, Y, phi[0], vmin=phi[0].min(), vmax=phi[0].max())
 text = ax.text(0.02, 1.02, "", transform=ax.transAxes)

examples/fixed_source/sphere_in_cube/input.py

@@ -0,0 +1,57 @@
+import numpy as np
+import mcdc
+
+# =============================================================================
+# Set model
+# =============================================================================
+# Homogeneous pure-fission sphere inside a pure-scattering cube
+
+# Set materials
+pure_f = mcdc.material(fission=np.array([1.0]), nu_p=np.array([1.2]))
+pure_s = mcdc.material(scatter=np.array([[1.0]]))
+
+# Set surfaces
+sx1 = mcdc.surface("plane-x", x=0.0, bc="vacuum")
+sx2 = mcdc.surface("plane-x", x=4.0, bc="vacuum")
+sy1 = mcdc.surface("plane-y", y=0.0, bc="vacuum")
+sy2 = mcdc.surface("plane-y", y=4.0, bc="vacuum")
+sz1 = mcdc.surface("plane-z", z=0.0, bc="vacuum")
+sz2 = mcdc.surface("plane-z", z=4.0, bc="vacuum")
+sphere = mcdc.surface("sphere", center=[2.0, 2.0, 2.0], radius=1.5)
+inside_sphere = -sphere
+inside_box = +sx1 & -sx2 & +sy1 & -sy2 & +sz1 & -sz2
+
+# Set cells
+# Source
+mcdc.cell(inside_box & ~inside_sphere, pure_s)
+
+# Sphere
+sphere_cell = mcdc.cell(inside_sphere, pure_f)
+
+# =============================================================================
+# Set source
+# =============================================================================
+# The source pulses in t=[0,5]
+
+mcdc.source(x=[0.0, 4.0], y=[0.0, 4.0], z=[0.0, 4.0], time=[0.0, 50.0], isotropic=True)
+
+# =============================================================================
+# Set tally, setting, and run mcdc
+# =============================================================================
+mcdc.tally.mesh_tally(
+    scores=["fission"],
+    x=np.linspace(0.0, 4.0, 2),
+    y=np.linspace(0.0, 4.0, 2),
+    z=np.linspace(0.0, 4.0, 2),
+    # t=np.linspace(0.0, 200.0, 2),
+)
+
+
+mcdc.tally.cell_tally(sphere_cell, scores=["fission"])
+
+# Setting
+mcdc.setting(N_particle=100)
+mcdc.implicit_capture()
+
+# Run
+mcdc.run()

examples/fixed_source/sphere_in_cube/process.py

@@ -0,0 +1,23 @@
+import h5py
+
+# Load results
+with h5py.File("output.h5", "r") as f:
+    # print(f["tallies"].keys())
+    print(f["input_deck"]["cell_tallies"].keys())
+
+    for i in range(len(f["input_deck"]["cell_tallies"])):
+        fission_score = f[f"tallies/cell_tally_{i}/fission"]
+
+        print(
+            f'for sphere {i+1}, mean = {fission_score["mean"][()]}, sdev = {fission_score["sdev"][()]}'
+        )
+
+        # print(fission_score["mean"][()])
+        # print(fission_score["sdev"][()])
+
+        # print(f"fission_score mean = {fission_score["mean"][()]}")
+        # print(f"fission_score mean = {fission_score["sdev"][()]}")
+
+    # cell = f["tallies/cell_tally_0/fission"]
+    # print(f'sphere1 mean = {cell["mean"][()]}')
+    # print(f'sphere2 sdev = {cell["sdev"][()]}')

mcdc/card.py

@@ -211,6 +211,8 @@ def __init__(self):
         self.fill_ID = None
         self.translation = np.array([0.0, 0.0, 0.0])
         self.surface_IDs = np.zeros(0, dtype=int)
+        self.N_tally = 0
+        self.tally_IDs = []
         self._region_RPN = []  # Reverse Polish Notation
 
     def set_region_RPN(self):
@@ -383,3 +385,12 @@ def __init__(self, surface_ID):
         # Set card data
         self.surface_ID = surface_ID
         self.N_bin = 1
+
+
+class CellTallyCard(TallyCard):
+    def __init__(self, cell_ID):
+        TallyCard.__init__(self, "Cell tally")
+
+        # Set card data
+        self.cell_ID = cell_ID
+        self.N_bin = 1