diff --git a/.gitignore b/.gitignore
index c3f1d70f..298e9ea2 100644
--- a/.gitignore
+++ b/.gitignore
@@ -37,3 +37,6 @@ tmp_*
 # test cache
 .pytest_cache
 pytestdebug.log
+
+# profiler
+*.prof
diff --git a/examples/c5g7/3d/TDX/input.py b/examples/c5g7/3d/TDX/input.py
index 54449490..593f9126 100644
--- a/examples/c5g7/3d/TDX/input.py
+++ b/examples/c5g7/3d/TDX/input.py
@@ -298,7 +298,7 @@ def set_mat(mat):
 mcdc.tally(scores=["fission"], t=t_grid)
 
 # Setting
-mcdc.setting(N_particle=1e4, active_bank_buff=1000)
+mcdc.setting(N_particle=1e2, active_bank_buff=1000)
 
 # Run
 mcdc.run()
diff --git a/examples/fixed_source/azurv1_pl_super/input.py b/examples/fixed_source/azurv1_pl_super/input.py
index b7ddd691..94e515d6 100644
--- a/examples/fixed_source/azurv1_pl_super/input.py
+++ b/examples/fixed_source/azurv1_pl_super/input.py
@@ -29,7 +29,7 @@
 # =============================================================================
 # Isotropic pulse at x=t=0
 
-mcdc.source(point=[0.0, 0.0, 0.0], isotropic=True)
+mcdc.source(point=[0.0, 0.0, 0.0], isotropic=True, time=[1e-10, 1e-10])
 
 # =============================================================================
 # Set tally, setting, and run mcdc
@@ -37,13 +37,13 @@
 
 # Tally: cell-average, cell-edge, and time-edge scalar fluxes
 mcdc.tally(
-    scores=["flux", "flux-x", "flux-t"],
+    scores=["flux"],
     x=np.linspace(-20.5, 20.5, 202),
     t=np.linspace(0.0, 20.0, 21),
 )
 
 # Setting
-mcdc.setting(N_particle=1e3)
+mcdc.setting(N_particle=1000)
 
 # Run
 mcdc.run()
diff --git a/examples/fixed_source/azurv1_pl_super/process.py b/examples/fixed_source/azurv1_pl_super/process.py
index 47eb1c74..85652134 100644
--- a/examples/fixed_source/azurv1_pl_super/process.py
+++ b/examples/fixed_source/azurv1_pl_super/process.py
@@ -3,46 +3,27 @@
 import matplotlib.animation as animation
 import h5py
 
-# =============================================================================
-# Reference solution (SS)
-# =============================================================================
-
-# Load grids
-with h5py.File("output.h5", "r") as f:
-    x = f["tally/grid/x"][:]
-    t = f["tally/grid/t"][:]
-
-dx = x[1] - x[0]
-x_mid = 0.5 * (x[:-1] + x[1:])
-dt = t[1:] - t[:-1]
-K = len(dt)
-J = len(x_mid)
 
+# Reference solution
 data = np.load("reference.npz")
-phi_x_ref = data["phi_x"]
-phi_t_ref = data["phi_t"]
 phi_ref = data["phi"]
 
-# =============================================================================
-# Animate results
-# =============================================================================
-
+# Get results
 with h5py.File("output.h5", "r") as f:
+    x = f["tally/grid/x"][:]
+    dx = x[1:] - x[:-1]
+    x_mid = 0.5 * (x[:-1] + x[1:])
+    t = f["tally/grid/t"][:]
+    dt = t[1:] - t[:-1]
+    K = len(t) - 1
+
     phi = f["tally/flux/mean"][:]
     phi_sd = f["tally/flux/sdev"][:]
-    phi_x = f["tally/flux-x/mean"][:]
-    phi_x_sd = f["tally/flux-x/sdev"][:]
-    phi_t = f["tally/flux-t/mean"][:]
-    phi_t_sd = f["tally/flux-t/sdev"][:]
-for k in range(K):
-    phi[k] /= dx * dt[k]
-    phi_sd[k] /= dx * dt[k]
-    phi_x[k] /= dt[k]
-    phi_x_sd[k] /= dt[k]
-    phi_t[k] /= dx
-    phi_t_sd[k] /= dx
-phi_t[K] /= dx
-phi_t_sd[K] /= dx
+
+    # Normalize
+    for k in range(K):
+        phi[k] /= dx * dt[k]
+        phi_sd[k] /= dx * dt[k]
 
 # Flux - average
 fig = plt.figure()
@@ -75,76 +56,3 @@ def animate(k):
 simulation = animation.FuncAnimation(fig, animate, frames=K)
 writervideo = animation.FFMpegWriter(fps=6)
 plt.show()
-
-# Flux - x
-fig = plt.figure()
-ax = plt.axes(
-    xlim=(-21.889999999999997, 21.89), ylim=(-0.042992644459595206, 0.9028455336514992)
-)
-ax.grid()
-ax.set_xlabel(r"$x$")
-ax.set_ylabel(r"Flux")
-ax.set_title(r"$\bar{\phi}_{k}(x)$")
-(line1,) = ax.plot([], [], "-b", label="MC")
-(line2,) = ax.plot([], [], "--r", label="Ref.")
-fb = ax.fill_between([], [], [], [], alpha=0.2, color="b")
-text = ax.text(0.02, 0.9, "", transform=ax.transAxes)
-ax.legend()
-
-
-def animate(k):
-    global fb
-    fb.remove()
-    line1.set_data(x, phi_x[k, :])
-    fb = ax.fill_between(
-        x,
-        phi_x[k, :] - phi_x_sd[k, :],
-        phi_x[k, :] + phi_x_sd[k, :],
-        alpha=0.2,
-        color="b",
-    )
-    line2.set_data(x, phi_x_ref[k, :])
-    text.set_text(r"$t \in [%.1f,%.1f]$ s" % (t[k], t[k + 1]))
-    return line1, line2, text
-
-
-simulation = animation.FuncAnimation(fig, animate, frames=K)
-writervideo = animation.FFMpegWriter(fps=6)
-plt.show()
-
-# Flux - t
-fig = plt.figure()
-ax = plt.axes(
-    xlim=(-21.889999999999997, 21.89), ylim=(-0.042992644459595206, 0.9028455336514992)
-)
-ax.grid()
-ax.set_xlabel(r"$x$")
-ax.set_ylabel(r"Flux")
-ax.set_title(r"$\bar{\phi}_{j}(t)$")
-(line1,) = ax.plot([], [], "-b", label="MC")
-(line2,) = ax.plot([], [], "--r", label="Ref.")
-fb = ax.fill_between([], [], [], [], alpha=0.2, color="b")
-text = ax.text(0.02, 0.9, "", transform=ax.transAxes)
-ax.legend()
-
-
-def animate(k):
-    global fb
-    fb.remove()
-    k += 1
-    line1.set_data(x_mid, phi_t[k, :])
-    fb = ax.fill_between(
-        x_mid,
-        phi_t[k, :] - phi_t_sd[k, :],
-        phi_t[k, :] + phi_t_sd[k, :],
-        alpha=0.2,
-        color="b",
-    )
-    line2.set_data(x_mid, phi_t_ref[k - 1, :])
-    text.set_text(r"$t = %.1f$ s" % (t[k]))
-    return line1, line2, text
-
-
-simulation = animation.FuncAnimation(fig, animate, frames=K)
-writervideo = animation.FFMpegWriter(fps=6)
-plt.show()
diff --git a/examples/fixed_source/azurv1_pl_super/reference.py b/examples/fixed_source/azurv1_pl_super/reference.py
index 120f5689..a2ad8f67 100644
--- a/examples/fixed_source/azurv1_pl_super/reference.py
+++ b/examples/fixed_source/azurv1_pl_super/reference.py
@@ -50,8 +50,6 @@ def phiX(x, t0, t1):
 
 
 phi_avg = np.zeros([K, J])
-phi_edge = np.zeros([K, J])
-phi_face = np.zeros([K, J + 1])
 
 for k in range(K):
     for j in range(J):
@@ -61,7 +59,6 @@ def phiX(x, t0, t1):
         t0 = t[k]
         t1 = t[k + 1]
         dt = t1 - t0
-        phi_edge[k, j] = quad(phi, x0, x1, args=(t1))[0] / dx
         phi_avg[k, j] = quad(phiX, x0, x1, args=(t0, t1))[0] / dx / dt
 
 for j in range(J + 1):
@@ -69,11 +66,7 @@ def phiX(x, t0, t1):
         t0 = t[k]
         t1 = t[k + 1]
         dt = t1 - t0
-        phi_face[k, j] = quad(phi_t, t0, t1, args=(x[j]))[0] / dt
 
-
-phi_edge = np.nan_to_num(phi_edge)
-phi_face = np.nan_to_num(phi_face)
 phi_avg = np.nan_to_num(phi_avg)
 
-np.savez("reference.npz", x=x, t=t, phi=phi_avg, phi_x=phi_face, phi_t=phi_edge)
+np.savez("reference.npz", x=x, t=t, phi=phi_avg)