From 0a523b5b90d73546aa70272fb63d9ef51f9d9509 Mon Sep 17 00:00:00 2001 From: Sam Pasmann Date: Thu, 9 Nov 2023 11:00:08 -0700 Subject: [PATCH] rename and add fission source tally functions --- mcdc/card.py | 1 - mcdc/input_.py | 2 -- mcdc/kernel.py | 42 ++++++++++++++++++++++++------------------ mcdc/loop.py | 4 +--- mcdc/type_.py | 3 +-- 5 files changed, 26 insertions(+), 26 deletions(-) diff --git a/mcdc/card.py b/mcdc/card.py index c05236bc..67d1d8ea 100644 --- a/mcdc/card.py +++ b/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]), diff --git a/mcdc/input_.py b/mcdc/input_.py index 6a05e7dd..b78cded1 100644 --- a/mcdc/input_.py +++ b/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: diff --git a/mcdc/kernel.py b/mcdc/kernel.py index 6a2e0721..259b8c0d 100644 --- a/mcdc/kernel.py +++ b/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"] diff --git a/mcdc/loop.py b/mcdc/loop.py index f817fa02..91393f36 100644 --- a/mcdc/loop.py +++ b/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 diff --git a/mcdc/type_.py b/mcdc/type_.py index e01842cc..803bdc7f 100644 --- a/mcdc/type_.py +++ b/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), ]