rename and add fission source tally functions

mcdc/card.py

@@ -96,7 +96,6 @@ def reset(self):
                 "maxitt": 5,
                 "N_dim": 6,
                 "seed": 12345,
-                "source_tilt": 0,
                 "scramble": False,
                 "fixed_source": np.ones([1, 1, 1, 1]),
                 "material_idx": np.ones([1, 1, 1, 1]),

mcdc/input_.py

@@ -1134,7 +1134,6 @@ def iQMC(
     N_dim=6,
     seed=12345,
     preconditioner_sweeps=5,
-    source_tilt=0,
     generator="halton",
     fixed_source_solver="source_iteration",
     eigenmode_solver="power_iteration",
@@ -1148,7 +1147,6 @@ def iQMC(
     card["iqmc"]["N_dim"] = N_dim
     card["iqmc"]["scramble"] = scramble
     card["iqmc"]["seed"] = seed
-    card["iqmc"]["source_tilt"] = source_tilt
 
     # Set mesh
     if g is not None:

mcdc/kernel.py

@@ -2324,11 +2324,7 @@ def iqmc_prepare_nuSigmaF(mcdc):
                     t = 0
                     mat_idx = mcdc["technique"]["iqmc"]["material_idx"][t, i, j, k]
                     material = mcdc["materials"][mat_idx]
-                    nu_f = material["nu_f"]
-                    SigmaF = material["fission"]
-                    mcdc["technique"]["iqmc"]["score"]["fission-source"][
-                        :, t, i, j, k
-                    ] = (nu_f * SigmaF * flux[:, t, i, j, k])
+                    mcdc["technique"]["iqmc"]["score"]["fission-source"] += iqmc_fission_source(flux[:, t, i, j, k], material)
 
 
 @njit
@@ -2358,10 +2354,10 @@ def iqmc_prepare_source(mcdc):
                 for k in range(Nz):
                     mat_idx = mcdc["technique"]["iqmc"]["material_idx"][t, i, j, k]
                     # we can vectorize the multigroup calculation here
-                    fission[:, t, i, j, k] = iqmc_fission_source(
+                    fission[:, t, i, j, k] = iqmc_effective_fission(
                         flux_fission[:, t, i, j, k], mat_idx, mcdc
                     )
-                    scatter[:, t, i, j, k] = iqmc_scattering_source(
+                    scatter[:, t, i, j, k] = iqmc_effective_scattering(
                         flux_scatter[:, t, i, j, k], mat_idx, mcdc
                     )
     mcdc["technique"]["iqmc"]["score"]["effective-scattering"] = scatter
@@ -2482,8 +2478,6 @@ def iqmc_score_tallies(P, distance, mcdc):
     material = mcdc["materials"][P["material_ID"]]
     w = P["iqmc"]["w"]
     SigmaT = material["total"]
-    SigmaF = material["fission"]
-    nu_f = material["nu_f"]
     mat_id = P["material_ID"]
 
     t, x, y, z, outside = mesh_get_index(P, mesh)
@@ -2506,18 +2500,18 @@ def iqmc_score_tallies(P, distance, mcdc):
     score_bin["flux"][:, t, x, y, z] += flux
 
     # Score effective source tallies
-    score_bin["effective-scattering"][:, t, x, y, z] += iqmc_scattering_source(
+    score_bin["effective-scattering"][:, t, x, y, z] += iqmc_effective_scattering(
         flux, mat_id, mcdc
     )
-    score_bin["effective-fission"][:, t, x, y, z] += iqmc_fission_source(
+    score_bin["effective-fission"][:, t, x, y, z] += iqmc_effective_fission(
         flux, mat_id, mcdc
     )
 
     if score_list["fission-source"]:
-        score_bin["fission-source"][:, t, x, y, z] += nu_f * SigmaF * flux
+        score_bin["fission-source"] += iqmc_fission_source(flux, material)
 
     if score_list["fission-power"]:
-        score_bin["fission-power"][:, t, x, y, z] += SigmaF * flux
+        score_bin["fission-power"][:, t, x, y, z] += iqmc_fission_power(flux, material)
 
     if score_list["tilt-t"]:
         t_mid = mesh["t"][t] + (dt * 0.5)
@@ -3023,7 +3017,20 @@ def iqmc_flux(SigmaT, w, distance, dV):
 
 
 @njit
-def iqmc_fission_source(phi, mat_id, mcdc):
+def iqmc_fission_source(phi, material):
+    SigmaF = material["fission"]
+    nu_f = material["nu_f"]
+    return np.sum(nu_f * SigmaF * phi)
+
+
+@njit
+def iqmc_fission_power(phi, material):
+    SigmaF = material["fission"]
+    return SigmaF * phi
+
+
+@njit
+def iqmc_effective_fission(phi, mat_id, mcdc):
     """
     Calculate the fission source for use with iQMC.
 
@@ -3058,7 +3065,7 @@ def iqmc_fission_source(phi, mat_id, mcdc):
 
 
 @njit
-def iqmc_scattering_source(phi, mat_id, mcdc):
+def iqmc_effective_scattering(phi, mat_id, mcdc):
     """
     Calculate the scattering source for use with iQMC.
 
@@ -3085,8 +3092,8 @@ def iqmc_scattering_source(phi, mat_id, mcdc):
 
 @njit
 def iqmc_effective_source(phi, mat_id, mcdc):
-    S = iqmc_scattering_source(phi, mat_id, mcdc)
-    F = iqmc_fission_source(phi, mat_id, mcdc)
+    S = iqmc_effective_scattering(phi, mat_id, mcdc)
+    F = iqmc_effective_fission(phi, mat_id, mcdc)
     return S + F
 
 
@@ -3274,7 +3281,6 @@ def FxV(V, mcdc):
     mcdc["bank_source"]["size"] = 0
 
     # QMC Sweep
-    # prepare_qmc_iqmc_fission_source(mcdc)
     mcdc["technique"]["iqmc"]["source"] = (
         mcdc["technique"]["iqmc"]["fixed_source"]
         + mcdc["technique"]["iqmc"]["score"]["effective-fission"]

mcdc/loop.py

@@ -542,8 +542,6 @@ def power_iteration(mcdc):
     fission_source_old = score_bin["fission-source"].copy()
 
     while not simulation_end:
-        # normalize source to effective fission source
-        # mcdc["technique"]["iqmc"]["source"] /= score_bin["effective-fission"].sum()
         # iterate over scattering source
         if solver == "source_iteration":
             mcdc["technique"]["iqmc"]["maxitt"] = 1
@@ -555,7 +553,7 @@ def power_iteration(mcdc):
         mcdc["technique"]["iqmc"]["itt"] = 0
 
         # update k_eff
-        mcdc["k_eff"] *= score_bin["fission-source"].sum() / fission_source_old.sum()
+        mcdc["k_eff"] *= score_bin["fission-source"] / fission_source_old
 
         # calculate diff in keff
         mcdc["technique"]["iqmc"]["res_outter"] = abs(mcdc["k_eff"] - k_old) / k_old

mcdc/type_.py

@@ -532,7 +532,7 @@ def make_type_technique(N_particle, G, card):
         ["tilt-yz", (Ng, Nt, Nx, Ny, Nz)],
         ["tilt-xyz", (Ng, Nt, Nx, Ny, Nz)],
         ["fission-power", (Ng, Nt, Nx, Ny, Nz)],  # SigmaF*phi
-        ["fission-source", (Ng, Nt, Nx, Ny, Nz)],  # nu*SigmaF*phi
+        ["fission-source", (1,)] # nu*SigmaF*phi
     ]
 
     if card["iQMC"]:
@@ -579,7 +579,6 @@ def make_type_technique(N_particle, G, card):
         ("krylov_restart", int64),
         ("preconditioner_sweeps", int64),
         ("sweep_counter", int64),
-        ("source_tilt", int64),
         ("w_min", float64),
     ]