Merge pull request CEMeNT-PSAAP#121 from ilhamv/purge_crossing_tally

remove crossing estimators

examples/fixed_source/kobayashi3-TD/input.py

@@ -68,7 +68,7 @@
 )
 
 # Setting
-mcdc.setting(N_particle=1e4)
+mcdc.setting(N_particle=1e6)
 mcdc.implicit_capture()
 
 # Run

examples/fixed_source/slab_absorbium/input.py

@@ -34,15 +34,15 @@
 # Set tally, setting, and run mcdc
 # =============================================================================
 
-# Tally: cell-average and cell-edge angular fluxes and currents
+# Tally: cell-average fluxes and currents
 mcdc.tally(
-    scores=["flux", "current", "flux-z", "current-z"],
+    scores=["flux", "current"],
     z=np.linspace(0.0, 6.0, 61),
     mu=np.linspace(-1.0, 1.0, 32 + 1),
 )
 
 # Setting
-mcdc.setting(N_particle=1e7)
+mcdc.setting(N_particle=1e3)
 
 # Run
 mcdc.run()

examples/fixed_source/slab_absorbium/process.py

@@ -17,28 +17,18 @@
 
     psi = f["tally/flux/mean"][:]
     psi_sd = f["tally/flux/sdev"][:]
-    psi_z = f["tally/flux-z/mean"][:]
-    psi_z_sd = f["tally/flux-z/sdev"][:]
     J = f["tally/current/mean"][:, 2]
     J_sd = f["tally/current/sdev"][:, 2]
-    J_z = f["tally/current-z/mean"][:, 2]
-    J_z_sd = f["tally/current-z/sdev"][:, 2]
 
 I = len(z) - 1
 N = len(mu) - 1
 
 # Scalar flux
 phi = np.zeros(I)
 phi_sd = np.zeros(I)
-phi_z = np.zeros(I + 1)
-phi_z_sd = np.zeros(I + 1)
 for i in range(I):
     phi[i] += np.sum(psi[i, :])
     phi_sd[i] += np.linalg.norm(psi_sd[i, :])
-    phi_z[i] += np.sum(psi_z[i, :])
-    phi_z_sd[i] += np.linalg.norm(psi_z_sd[i, :])
-phi_z[I] += np.sum(psi_z[I, :])
-phi_z_sd[I] += np.linalg.norm(psi_z_sd[I, :])
 
 # Normalize
 phi /= dz
@@ -50,7 +40,7 @@
     psi_sd[:, n] = psi_sd[:, n] / dz / dmu[n]
 
 # Reference solution
-phi_ref, phi_z_ref, J_ref, J_z_ref, psi_ref, psi_z_ref = reference(z, mu)
+phi_ref, J_ref, psi_ref = reference(z, mu)
 
 # Flux - spatial average
 plt.plot(z_mid, phi, "-b", label="MC")
@@ -64,18 +54,6 @@
 plt.title(r"$\bar{\phi}_i$")
 plt.show()
 
-# Flux - spatial grid
-plt.plot(z, phi_z, "-b", label="MC")
-plt.fill_between(z, phi_z - phi_z_sd, phi_z + phi_z_sd, alpha=0.2, color="b")
-plt.plot(z, phi_z_ref, "--r", label="Ref.")
-plt.xlabel(r"$z$, cm")
-plt.ylabel("Flux")
-plt.ylim([0.06, 0.16])
-plt.grid()
-plt.legend()
-plt.title(r"$\phi(z)$")
-plt.show()
-
 # Current - spatial average
 plt.plot(z_mid, J, "-b", label="MC")
 plt.fill_between(z_mid, J - J_sd, J + J_sd, alpha=0.2, color="b")
@@ -88,18 +66,6 @@
 plt.title(r"$\bar{J}_i$")
 plt.show()
 
-# Current - spatial grid
-plt.plot(z, J_z, "-b", label="MC")
-plt.fill_between(z, J_z - J_z_sd, J_z + J_z_sd, alpha=0.2, color="b")
-plt.plot(z, J_z_ref, "--r", label="Ref.")
-plt.xlabel(r"$z$, cm")
-plt.ylabel("Current")
-plt.ylim([-0.03, 0.045])
-plt.grid()
-plt.legend()
-plt.title(r"$J(z)$")
-plt.show()
-
 # Angular flux - spatial average
 vmin = min(np.min(psi_ref), np.min(psi))
 vmax = max(np.max(psi_ref), np.max(psi))

examples/fixed_source/slab_absorbium/reference.py

@@ -101,45 +101,30 @@ def psi3_(x, mu0, mu1):
     phi = np.zeros(I)
     psi = np.zeros((I, N))
     J = np.zeros(I)
-    phi_x = np.zeros(I + 1)
-    psi_x = np.zeros((I + 1, N))
-    J_x = np.zeros(I + 1)
 
     for i in range(int(I / 3)):
         phi[i] = quad(phi1, x[i], x[i + 1])[0] / dx[i]
         J[i] = quad(J1, x[i], x[i + 1])[0] / dx[i]
-        phi_x[i] = phi1(x[i])
-        J_x[i] = J1(x[i])
         for n in range(N):
             mu0 = mu[n]
             mu1 = mu[n + 1]
             psi[i, n] = quad(psi1_, x[i], x[i + 1], args=(mu0, mu1))[0] / dx[i] / dmu[n]
-            psi_x[i, n] = psi1_(x[i], mu0, mu1) / dmu[n]
     for i in range(int(I / 3), int(2 * I / 3)):
         phi[i] = quad(phi2, x[i], x[i + 1])[0] / dx[i]
         J[i] = quad(J2, x[i], x[i + 1])[0] / dx[i]
-        phi_x[i] = phi2(x[i])
-        J_x[i] = J2(x[i])
         for n in range(N):
             mu0 = mu[n]
             mu1 = mu[n + 1]
             psi[i, n] = quad(psi2_, x[i], x[i + 1], args=(mu0, mu1))[0] / dx[i] / dmu[n]
-            psi_x[i, n] = psi2_(x[i], mu0, mu1) / dmu[n]
     for i in range(int(2 * I / 3), I):
         phi[i] = quad(phi3, x[i], x[i + 1])[0] / dx[i]
         J[i] = quad(J3, x[i], x[i + 1])[0] / dx[i]
-        phi_x[i] = phi3(x[i])
-        J_x[i] = J3(x[i])
         for n in range(N):
             mu0 = mu[n]
             mu1 = mu[n + 1]
             psi[i, n] = quad(psi3_, x[i], x[i + 1], args=(mu0, mu1))[0] / dx[i] / dmu[n]
-            psi_x[i, n] = psi3_(x[i], mu0, mu1) / dmu[n]
-    phi_x[I] = phi3(x[I])
-    J_x[I] = J3(x[I])
     for n in range(N):
         mu0 = mu[n]
         mu1 = mu[n + 1]
-        psi_x[I, n] = psi3_(x[I], mu0, mu1) / dmu[n]
 
-    return phi, phi_x, J, J_x, psi, psi_x
+    return phi, J, psi

mcdc/card.py

@@ -24,42 +24,13 @@ def reset(self):
 
         self.tally = {
             "tag": "Tally",
-            "tracklength": False,
+            "tracklength": True,
             "flux": False,
             "density": False,
             "fission": False,
             "total": False,
             "current": False,
             "eddington": False,
-            "crossing": False,
-            "crossing_x": False,
-            "flux_x": False,
-            "density_x": False,
-            "fission_x": False,
-            "total_x": False,
-            "current_x": False,
-            "eddington_x": False,
-            "crossing_y": False,
-            "flux_y": False,
-            "density_y": False,
-            "fission_y": False,
-            "total_y": False,
-            "current_y": False,
-            "eddington_y": False,
-            "crossing_z": False,
-            "flux_z": False,
-            "density_z": False,
-            "fission_z": False,
-            "total_z": False,
-            "current_z": False,
-            "eddington_z": False,
-            "crossing_t": False,
-            "flux_t": False,
-            "density_t": False,
-            "fission_t": False,
-            "total_t": False,
-            "current_t": False,
-            "eddington_t": False,
             "mesh": make_card_mesh(),
         }
 

mcdc/constant.py

@@ -15,12 +15,6 @@
 EVENT_LATTICE = 1 << 8
 EVENT_SURFACE_MOVE = 1 << 9
 
-# Mesh crossing flags
-MESH_X = 0
-MESH_Y = 1
-MESH_Z = 2
-MESH_T = 3
-
 # Gyration raius type
 GYRATION_RADIUS_ALL = 0
 GYRATION_RADIUS_INFINITE_X = 1

mcdc/input_.py

@@ -856,38 +856,13 @@ def tally(
         found = False
         for score_name in type_.score_list:
             if s.replace("-", "_") == score_name:
+                card["tracklength"] = True
                 card[score_name] = True
                 found = True
                 break
         if not found:
             print_error("Unknown tally score %s" % s)
 
-    # Set estimator flags
-    for score_name in type_.score_tl_list:
-        if card[score_name]:
-            card["tracklength"] = True
-            break
-    for score_name in type_.score_x_list:
-        if card[score_name]:
-            card["crossing"] = True
-            card["crossing_x"] = True
-            break
-    for score_name in type_.score_y_list:
-        if card[score_name]:
-            card["crossing"] = True
-            card["crossing_y"] = True
-            break
-    for score_name in type_.score_z_list:
-        if card[score_name]:
-            card["crossing"] = True
-            card["crossing_z"] = True
-            break
-    for score_name in type_.score_t_list:
-        if card[score_name]:
-            card["crossing"] = True
-            card["crossing_t"] = True
-            break
-
     return card