diff --git a/mcdc/__init__.py b/mcdc/__init__.py index c902b276..c36bea2c 100644 --- a/mcdc/__init__.py +++ b/mcdc/__init__.py @@ -18,7 +18,6 @@ iQMC, weight_roulette, IC_generator, - dsm, uq, reset, domain_decomposition, diff --git a/mcdc/card.py b/mcdc/card.py index fa53a1fa..3a989e55 100644 --- a/mcdc/card.py +++ b/mcdc/card.py @@ -54,9 +54,6 @@ def __init__(self, G=1, J=0): self.chi_s = np.zeros([G, G]) self.chi_p = np.zeros([G, G]) self.chi_d = np.zeros([J, G]) - self.sensitivity = False - self.sensitivity_ID = 0 - self.dsm_Np = 1.0 self.uq = False self.flags = [] self.distribution = "" @@ -85,7 +82,6 @@ def __init__(self, N_nuclide, G=1, J=0): self.nu_f = np.zeros(G) self.chi_s = np.zeros([G, G]) self.chi_p = np.zeros([G, G]) - self.sensitivity = False self.uq = False self.flags = [] self.distribution = "" @@ -167,9 +163,6 @@ def __init__(self): self.nx = 0.0 self.ny = 0.0 self.nz = 0.0 - self.sensitivity = False - self.sensitivity_ID = 0 - self.dsm_Np = 1.0 def _create_halfspace(self, positive): region = RegionCard("halfspace") diff --git a/mcdc/global_.py b/mcdc/global_.py index 7de89fbb..8df838f7 100644 --- a/mcdc/global_.py +++ b/mcdc/global_.py @@ -78,8 +78,6 @@ def reset(self): # Below are parameters not copied to mcdc.setting "bank_active_buff": 100, "bank_census_buff": 1.0, - # TODO: Move to technique - "N_sensitivity": 0, } self.technique = { @@ -158,7 +156,6 @@ def reset(self): "IC_precursor_density_max": 0.0, "IC_cycle_stretch": 1.0, "branchless_collision": False, - "dsm_order": 1, "uq": False, } diff --git a/mcdc/input_.py b/mcdc/input_.py index 02c14c77..b40e92f3 100644 --- a/mcdc/input_.py +++ b/mcdc/input_.py @@ -52,8 +52,6 @@ def nuclide( chi_d=None, speed=None, decay=None, - sensitivity=False, - dsm_Np=1.0, ): """ Create a nuclide @@ -80,11 +78,6 @@ def nuclide( Energy group speed [cm/s]. decay : numpy.ndarray (1D), optional Precursor group decay constant [/s]. - sensitivity : bool, optional - Set to `True` to calculate sensitivities to the nuclide. - dsm_Np : float - Average number of derivative particles produced at each - sensitivity nuclide collision. Returns ------- @@ -208,20 +201,6 @@ def nuclide( if np.sum(card.chi_d[dg, :]) > 0.0: card.chi_d[dg, :] /= np.sum(card.chi_d[dg, :]) - # Sensitivity setup - if sensitivity: - # Set flag - card.sensitivity = True - global_.input_deck.technique["sensitivity"] = True - global_.input_deck.technique["weighted_emission"] = False - - # Set ID - global_.input_deck.setting["N_sensitivity"] += 1 - card.sensitivity_ID = global_.input_deck.setting["N_sensitivity"] - - # Set dsm_Np - card.dsm_Np = dsm_Np - # Add to deck global_.input_deck.nuclides.append(card) @@ -240,8 +219,6 @@ def material( chi_d=None, speed=None, decay=None, - sensitivity=False, - dsm_Np=1.0, ): """ Create a material @@ -273,12 +250,6 @@ def material( Energy group speed [cm/s]. decay : numpy.ndarray (1D), optional Precursor group decay constant [/s]. - sensitivity : bool, optional - Set to `True` to calculate sensitivities to the material - (only relevant for single-nuclide material). - dsm_Np : float - Average number of derivative particles produced at each - sensitivity material collision (only relevant for single_nuclide material). Returns ------- @@ -303,8 +274,6 @@ def material( chi_d, speed, decay, - sensitivity, - dsm_Np, ) nuclides = [[card_nuclide, 1.0]] @@ -379,7 +348,7 @@ def material( # Set ID card.ID = len(global_.input_deck.materials) - # Calculate basic XS and determine sensitivity flag + # Calculate basic XS for i in range(N_nuclide): nuc = nuclides[i][0] density = nuclides[i][1] @@ -390,8 +359,6 @@ def material( card.scatter += nuc.scatter * density card.fission += nuc.fission * density card.total += nuc.total * density - card.sensitivity += nuc.sensitivity * density - card.sensitivity = bool(card.sensitivity) # Calculate effective speed # Current approach: weighted by nuclide macroscopic total cross section @@ -450,7 +417,7 @@ def material( return card -def surface(type_, bc="interface", sensitivity=False, dsm_Np=1.0, **kw): +def surface(type_, bc="interface", **kw): """ Create a surface to define the region of a cell. @@ -461,11 +428,6 @@ def surface(type_, bc="interface", sensitivity=False, dsm_Np=1.0, **kw): Surface type. bc : {"interface", "vacuum", "reflective"} Surface boundary condition. - sensitivity : bool, optional - Set to `True` to calculate sensitivities to the surface position. - dsm_Np : int - Average number of derivative particles produced at each - sensitivity surface crossing. Other Parameters ---------------- @@ -538,20 +500,6 @@ def surface(type_, bc="interface", sensitivity=False, dsm_Np=1.0, **kw): ) card.boundary_type = bc - # Sensitivity - if sensitivity: - # Set flag - card.sensitivity = True - global_.input_deck.technique["sensitivity"] = True - global_.input_deck.technique["weighted_emission"] = False - - # Set ID - global_.input_deck.setting["N_sensitivity"] += 1 - card.sensitivity_ID = global_.input_deck.setting["N_sensitivity"] - - # Set dsm_Np - card.dsm_Np = dsm_Np - # ========================================================================== # Surface attributes # ========================================================================== @@ -1776,22 +1724,6 @@ def IC_generator( card_setting["N_active"] = N_cycle -def dsm(order=1): - """ - Direct sensitivity method - - Parameters - ---------- - order : int, optional - order of the sensitivity to probe, by default 1 - """ - - card = global_.input_deck.technique - if order > 2: - print_error("DSM currently only supports up to second-order sensitivities") - card["dsm_order"] = order - - def uq(**kw): """ Activate uncertainty quantification. @@ -1950,7 +1882,6 @@ def make_particle_bank(size): ("g", np.uint64), ("E", np.float64), ("w", np.float64), - ("sensitivity_ID", np.int64), ("rng_seed", np.uint64), ] iqmc_struct = [("w", np.float64, (1,))] diff --git a/mcdc/kernel.py b/mcdc/kernel.py index a02f321f..8992beb0 100644 --- a/mcdc/kernel.py +++ b/mcdc/kernel.py @@ -654,7 +654,6 @@ def source_particle_dd(seed, mcdc): P["uz"] = uz P["g"] = g P["w"] = 1 - P["sensitivity_ID"] = 0 return P @@ -909,8 +908,6 @@ def source_particle(seed, mcdc): P["E"] = E P["w"] = 1.0 - P["sensitivity_ID"] = 0 - return P @@ -987,7 +984,6 @@ def get_particle(P, bank, mcdc): P["iqmc"]["w"] = P_rec["iqmc"]["w"] P["alive"] = True - P["sensitivity_ID"] = P_rec["sensitivity_ID"] # Set default IDs and event P["material_ID"] = -1 @@ -1644,7 +1640,6 @@ def copy_recordlike(P_new, P): P_new["E"] = P["E"] P_new["w"] = P["w"] P_new["rng_seed"] = P["rng_seed"] - P_new["sensitivity_ID"] = P["sensitivity_ID"] P_new["iqmc"]["w"] = P["iqmc"]["w"] copy_track_data(P_new, P) @@ -1687,7 +1682,6 @@ def copy_particle(P_new, P): P_new["surface_ID"] = P["surface_ID"] P_new["translation"] = P["translation"] P_new["event"] = P["event"] - P_new["sensitivity_ID"] = P["sensitivity_ID"] P_new["rng_seed"] = P["rng_seed"] P_new["iqmc"]["w"] = P["iqmc"]["w"] copy_track_data(P_new, P) @@ -2175,7 +2169,6 @@ def score_tracklength(P, distance, mcdc): material = mcdc["materials"][P["material_ID"]] # Get indices - s = P["sensitivity_ID"] t, x, y, z, outside = mesh_get_index(P, tally["mesh"]) mu, azi = mesh_get_angular_index(P, tally["mesh"]) g, outside_energy = mesh_get_energy_index(P, tally["mesh"], mcdc) @@ -2187,20 +2180,20 @@ def score_tracklength(P, distance, mcdc): # Score flux = distance * P["w"] if tally["flux"]: - score_flux(s, g, t, x, y, z, mu, azi, flux, tally["score"]["flux"]) + score_flux(g, t, x, y, z, mu, azi, flux, tally["score"]["flux"]) if tally["density"]: flux /= get_particle_speed(P, mcdc) - score_flux(s, g, t, x, y, z, mu, azi, flux, tally["score"]["density"]) + score_flux(g, t, x, y, z, mu, azi, flux, tally["score"]["density"]) if tally["fission"]: flux *= get_MacroXS(XS_FISSION, material, P, mcdc) - score_flux(s, g, t, x, y, z, mu, azi, flux, tally["score"]["fission"]) + score_flux(g, t, x, y, z, mu, azi, flux, tally["score"]["fission"]) if tally["total"]: flux *= get_MacroXS(XS_TOTAL, material, P, mcdc) - score_flux(s, g, t, x, y, z, mu, azi, flux, tally["score"]["total"]) + score_flux(g, t, x, y, z, mu, azi, flux, tally["score"]["total"]) if tally["current"]: - score_current(s, g, t, x, y, z, flux, P, tally["score"]["current"]) + score_current(g, t, x, y, z, flux, P, tally["score"]["current"]) if tally["eddington"]: - score_eddington(s, g, t, x, y, z, flux, P, tally["score"]["eddington"]) + score_eddington(g, t, x, y, z, flux, P, tally["score"]["eddington"]) @njit @@ -2208,7 +2201,6 @@ def score_exit(P, x, mcdc): tally = mcdc["tally"] material = mcdc["materials"][P["material_ID"]] - s = P["sensitivity_ID"] mu, azi = mesh_get_angular_index(P, tally["mesh"]) g, outside_energy = mesh_get_energy_index(P, tally["mesh"], mcdc) @@ -2218,27 +2210,27 @@ def score_exit(P, x, mcdc): # Score flux = P["w"] / abs(P["ux"]) - score_flux(s, g, 0, x, 0, 0, mu, azi, flux, tally["score"]["exit"]) + score_flux(g, 0, x, 0, 0, mu, azi, flux, tally["score"]["exit"]) @njit -def score_flux(s, g, t, x, y, z, mu, azi, flux, score): +def score_flux(g, t, x, y, z, mu, azi, flux, score): # score["bin"][s, g, t, x, y, z, mu, azi] += flux - adapt.global_add(score["bin"], (s, g, t, x, y, z, mu, azi), flux) + adapt.global_add(score["bin"], (g, t, x, y, z, mu, azi), flux) @njit -def score_current(s, g, t, x, y, z, flux, P, score): +def score_current(g, t, x, y, z, flux, P, score): # score["bin"][s, g, t, x, y, z, 0] += flux * P["ux"] # score["bin"][s, g, t, x, y, z, 1] += flux * P["uy"] # score["bin"][s, g, t, x, y, z, 2] += flux * P["uz"] - adapt.global_add(score["bin"], (s, g, t, x, y, z, 0), flux * P["ux"]) - adapt.global_add(score["bin"], (s, g, t, x, y, z, 1), flux * P["uy"]) - adapt.global_add(score["bin"], (s, g, t, x, y, z, 2), flux * P["uz"]) + adapt.global_add(score["bin"], (g, t, x, y, z, 0), flux * P["ux"]) + adapt.global_add(score["bin"], (g, t, x, y, z, 1), flux * P["uy"]) + adapt.global_add(score["bin"], (g, t, x, y, z, 2), flux * P["uz"]) @njit -def score_eddington(s, g, t, x, y, z, flux, P, score): +def score_eddington(g, t, x, y, z, flux, P, score): ux = P["ux"] uy = P["uy"] uz = P["uz"] @@ -2248,12 +2240,12 @@ def score_eddington(s, g, t, x, y, z, flux, P, score): # score["bin"][s, g, t, x, y, z, 3] += flux * uy * uy # score["bin"][s, g, t, x, y, z, 4] += flux * uy * uz # score["bin"][s, g, t, x, y, z, 5] += flux * uz * uz - adapt.global_add(score["bin"], (s, g, t, x, y, z, 0), flux * ux * ux) - adapt.global_add(score["bin"], (s, g, t, x, y, z, 1), flux * ux * uy) - adapt.global_add(score["bin"], (s, g, t, x, y, z, 2), flux * ux * uz) - adapt.global_add(score["bin"], (s, g, t, x, y, z, 3), flux * uy * uy) - adapt.global_add(score["bin"], (s, g, t, x, y, z, 4), flux * uy * uz) - adapt.global_add(score["bin"], (s, g, t, x, y, z, 5), flux * uz * uz) + adapt.global_add(score["bin"], (g, t, x, y, z, 0), flux * ux * ux) + adapt.global_add(score["bin"], (g, t, x, y, z, 1), flux * ux * uy) + adapt.global_add(score["bin"], (g, t, x, y, z, 2), flux * ux * uz) + adapt.global_add(score["bin"], (g, t, x, y, z, 3), flux * uy * uy) + adapt.global_add(score["bin"], (g, t, x, y, z, 4), flux * uy * uz) + adapt.global_add(score["bin"], (g, t, x, y, z, 5), flux * uz * uz) @njit @@ -2800,9 +2792,6 @@ def surface_crossing(P, prog): # Small shift to ensure crossing surface_shift(P, surface, trans, mcdc) - # Record old material for sensitivity quantification - material_ID_old = P["material_ID"] - # Tally particle exit if mcdc["tally"]["exit"] and not P["alive"]: # Reflectance if P["surface_ID"] == 0, else transmittance @@ -2818,17 +2807,6 @@ def surface_crossing(P, prog): trans = trans_struct["values"] P["cell_ID"] = get_particle_cell(P, UNIVERSE_ROOT, trans, mcdc) - # Sensitivity quantification for surface? - if surface["sensitivity"] and ( - P["sensitivity_ID"] == 0 - or mcdc["technique"]["dsm_order"] == 2 - and P["sensitivity_ID"] <= mcdc["setting"]["N_sensitivity"] - ): - material_ID_new = get_particle_material(P, mcdc) - if material_ID_old != material_ID_new: - # Sample derivative source particles - sensitivity_surface(P, surface, material_ID_old, material_ID_new, prog) - # ============================================================================= # Collision @@ -2940,7 +2918,6 @@ def sample_phasespace_scattering(P, material, P_new, mcdc): P_new["y"] = P["y"] P_new["z"] = P["z"] P_new["t"] = P["t"] - P_new["sensitivity_ID"] = P["sensitivity_ID"] if mcdc["setting"]["mode_MG"]: scattering_MG(P, material, P_new) @@ -2955,7 +2932,6 @@ def sample_phasespace_scattering_nuclide(P, nuclide, P_new, mcdc): P_new["y"] = P["y"] P_new["z"] = P["z"] P_new["t"] = P["t"] - P_new["sensitivity_ID"] = P["sensitivity_ID"] scattering_MG(P, nuclide, P_new) @@ -3230,7 +3206,6 @@ def sample_phasespace_fission(P, material, P_new, mcdc): P_new["y"] = P["y"] P_new["z"] = P["z"] P_new["t"] = P["t"] - P_new["sensitivity_ID"] = P["sensitivity_ID"] # Sample isotropic direction P_new["ux"], P_new["uy"], P_new["uz"] = sample_isotropic_direction(P_new) @@ -3291,7 +3266,6 @@ def sample_phasespace_fission_nuclide(P, nuclide, P_new, mcdc): P_new["y"] = P["y"] P_new["z"] = P["z"] P_new["t"] = P["t"] - P_new["sensitivity_ID"] = P["sensitivity_ID"] # Sample isotropic direction P_new["ux"], P_new["uy"], P_new["uz"] = sample_isotropic_direction(P_new) @@ -3516,346 +3490,6 @@ def weight_roulette(P, mcdc): P["alive"] = False -# ============================================================================== -# Sensitivity quantification (Derivative Source Method) -# ============================================================================= - - -@njit -def sensitivity_surface(P, surface, material_ID_old, material_ID_new, prog): - - mcdc = adapt.device(prog) - - # Sample number of derivative sources - xi = surface["dsm_Np"] - if xi != 1.0: - Np = int(math.floor(xi + rng(P))) - else: - Np = 1 - - # Terminate and put the current particle into the secondary bank - P["alive"] = False - adapt.add_active(copy_record(P), prog) - - # Get sensitivity ID - ID = surface["sensitivity_ID"] - if mcdc["technique"]["dsm_order"] == 2: - ID1 = min(P["sensitivity_ID"], ID) - ID2 = max(P["sensitivity_ID"], ID) - ID = get_DSM_ID(ID1, ID2, mcdc["setting"]["N_sensitivity"]) - - # Get materials - material_old = mcdc["materials"][material_ID_old] - material_new = mcdc["materials"][material_ID_new] - - # Determine the plus and minus components and then their weight signs - trans = P["translation"] - sign_origin = surface_normal_component(P, surface, trans) - if sign_origin > 0.0: - # New is +, old is - - sign_new = -1.0 - sign_old = 1.0 - else: - sign_new = 1.0 - sign_old = -1.0 - - # Get XS - g = P["g"] - SigmaT_old = material_old["total"][g] - SigmaT_new = material_new["total"][g] - SigmaS_old = material_old["scatter"][g] - SigmaS_new = material_new["scatter"][g] - SigmaF_old = material_old["fission"][g] - SigmaF_new = material_new["fission"][g] - nu_s_old = material_old["nu_s"][g] - nu_s_new = material_new["nu_s"][g] - nu_old = material_old["nu_f"][g] - nu_new = material_new["nu_f"][g] - nuSigmaS_old = nu_s_old * SigmaS_old - nuSigmaS_new = nu_s_new * SigmaS_new - nuSigmaF_old = nu_old * SigmaF_old - nuSigmaF_new = nu_new * SigmaF_new - - # Get source type probabilities - delta = -(SigmaT_old * sign_old + SigmaT_new * sign_new) - scatter = nuSigmaS_old * sign_old + nuSigmaS_new * sign_new - fission = nuSigmaF_old * sign_old + nuSigmaF_new * sign_new - p_delta = abs(delta) - p_scatter = abs(scatter) - p_fission = abs(fission) - p_total = p_delta + p_scatter + p_fission - - # Get inducing flux - # Apply constant flux approximation for tangent direction - # [Dupree 2002, Eq. (7.39)] - mu = abs(sign_origin) - epsilon = 0.01 - if mu < epsilon: - mu = epsilon / 2 - flux = P["w"] / mu - - # Base weight - w_hat = p_total * flux / xi - - # Sample the derivative sources - for n in range(Np): - # Create new particle - P_new = split_particle(P) - - # Sample source type - xi = rng(P) * p_total - tot = p_delta - if tot > xi: - # Delta source - sign_delta = delta / p_delta - P_new["w"] = w_hat * sign_delta - else: - tot += p_scatter - if tot > xi: - # Scattering source - total_scatter = nuSigmaS_old + nuSigmaS_new - w_s = w_hat * total_scatter / p_scatter - - # Sample if it is from + or - component - if nuSigmaS_old > rng(P) * total_scatter: - sample_phasespace_scattering(P, material_old, P_new, mcdc) - P_new["w"] = w_s * sign_old - else: - sample_phasespace_scattering(P, material_new, P_new, mcdc) - P_new["w"] = w_s * sign_new - else: - # Fission source - total_fission = nuSigmaF_old + nuSigmaF_new - w_f = w_hat * total_fission / p_fission - - # Sample if it is from + or - component - if nuSigmaF_old > rng(P) * total_fission: - sample_phasespace_fission(P, material_old, P_new, mcdc) - P_new["w"] = w_f * sign_old - else: - sample_phasespace_fission(P, material_new, P_new, mcdc) - P_new["w"] = w_f * sign_new - - # Assign sensitivity_ID - P_new["sensitivity_ID"] = ID - - # Shift back if needed to ensure crossing - sign = surface_normal_component(P_new, surface, trans) - if sign_origin * sign > 0.0: - # Get surface normal - nx, ny, nz = surface_normal(P_new, surface, trans) - - # The shift - if sign > 0.0: - P_new["x"] -= nx * 2 * SHIFT - P_new["y"] -= ny * 2 * SHIFT - P_new["z"] -= nz * 2 * SHIFT - else: - P_new["x"] += nx * 2 * SHIFT - P_new["y"] += ny * 2 * SHIFT - P_new["z"] += nz * 2 * SHIFT - - # Put the current particle into the secondary bank - adapt.add_active(P_new, prog) - - # Sample potential second-order sensitivity particles? - if mcdc["technique"]["dsm_order"] < 2 or P["sensitivity_ID"] > 0: - return - - # Get total probability - p_total = 0.0 - for material in [material_new, material_old]: - if material["sensitivity"]: - N_nuclide = material["N_nuclide"] - for i in range(N_nuclide): - nuclide = mcdc["nuclides"][material["nuclide_IDs"][i]] - if nuclide["sensitivity"]: - sigmaT = nuclide["total"][g] - sigmaS = nuclide["scatter"][g] - sigmaF = nuclide["fission"][g] - nu_s = nuclide["nu_s"][g] - nu = nuclide["nu_f"][g] - nusigmaS = nu_s * sigmaS - nusigmaF = nu * sigmaF - total = sigmaT + nusigmaS + nusigmaF - p_total += total - - # Base weight - w = p_total * flux / surface["dsm_Np"] - - # Sample source - for n in range(Np): - source_obtained = False - - # Create new particle - P_new = split_particle(P) - - # Sample term - xi = rng(P_new) * p_total - tot = 0.0 - - for idx in range(2): - - if idx == 0: - material_ID = material_ID_old - sign = sign_old - else: - material_ID = material_ID_new - sign = sign_new - - material = mcdc["materials"][material_ID] - if material["sensitivity"]: - N_nuclide = material["N_nuclide"] - for i in range(N_nuclide): - nuclide = mcdc["nuclides"][material["nuclide_IDs"][i]] - if nuclide["sensitivity"]: - # Source ID - ID1 = min(nuclide["sensitivity_ID"], surface["sensitivity_ID"]) - ID2 = max(nuclide["sensitivity_ID"], surface["sensitivity_ID"]) - ID_source = get_DSM_ID( - ID1, ID2, mcdc["setting"]["N_sensitivity"] - ) - - sigmaT = nuclide["total"][g] - sigmaS = nuclide["scatter"][g] - sigmaF = nuclide["fission"][g] - nu_s = nuclide["nu_s"][g] - nu = nuclide["nu_f"][g] - nusigmaS = nu_s * sigmaS - nusigmaF = nu * sigmaF - - tot += sigmaT - if tot > xi: - # Delta source - P_new["w"] = -w * sign - P_new["sensitivity_ID"] = ID_source - adapt.add_active(P_new, prog) - source_obtained = True - else: - P_new["w"] = w * sign - - tot += nusigmaS - if tot > xi: - # Scattering source - sample_phasespace_scattering_nuclide( - P, nuclide, P_new, mcdc - ) - P_new["sensitivity_ID"] = ID_source - adapt.add_active(P_new, prog) - source_obtained = True - else: - tot += nusigmaF - if tot > xi: - # Fission source - sample_phasespace_fission_nuclide( - P, nuclide, P_new, mcdc - ) - P_new["sensitivity_ID"] = ID_source - adapt.add_active(P_new, prog) - source_obtained = True - if source_obtained: - break - if source_obtained: - break - - -@njit -def sensitivity_material(P, prog): - - mcdc = adapt.device(prog) - - # The incident particle is already terminated - - # Get material - material = mcdc["materials"][P["material_ID"]] - - # Check if sensitivity nuclide is sampled - g = P["g"] - SigmaT = material["total"][g] - N_nuclide = material["N_nuclide"] - if N_nuclide == 1: - nuclide = mcdc["nuclides"][material["nuclide_IDs"][0]] - else: - xi = rng(P) * SigmaT - tot = 0.0 - for i in range(N_nuclide): - nuclide = mcdc["nuclides"][material["nuclide_IDs"][i]] - density = material["nuclide_densities"][i] - tot += density * nuclide["total"][g] - if xi < tot: - break - if not nuclide["sensitivity"]: - return - - # Sample number of derivative sources - xi = nuclide["dsm_Np"] - if xi != 1.0: - Np = int(math.floor(xi + rng(P))) - else: - Np = 1 - - # Get sensitivity ID - ID = nuclide["sensitivity_ID"] - double = False - if mcdc["technique"]["dsm_order"] == 2: - ID1 = min(P["sensitivity_ID"], ID) - ID2 = max(P["sensitivity_ID"], ID) - ID = get_DSM_ID(ID1, ID2, mcdc["setting"]["N_sensitivity"]) - if ID1 == ID2: - double = True - - # Undo implicit capture - if mcdc["technique"]["implicit_capture"]: - SigmaC = material["capture"][g] - P["w"] *= SigmaT / (SigmaT - SigmaC) - - # Get XS - g = P["g"] - sigmaT = nuclide["total"][g] - sigmaS = nuclide["scatter"][g] - sigmaF = nuclide["fission"][g] - nu_s = nuclide["nu_s"][g] - nu = nuclide["nu_f"][g] - nusigmaS = nu_s * sigmaS - nusigmaF = nu * sigmaF - - # Base weight - total = sigmaT + nusigmaS + nusigmaF - w = total * P["w"] / sigmaT / xi - - # Double if it's self-second-order - if double: - w *= 2 - - # Sample the derivative sources - for n in range(Np): - # Create new particle - P_new = split_particle(P) - - # Sample source type - xi = rng(P_new) * total - tot = sigmaT - if tot > xi: - # Delta source - P_new["w"] = -w - else: - P_new["w"] = w - - tot += nusigmaS - if tot > xi: - # Scattering source - sample_phasespace_scattering_nuclide(P, nuclide, P_new, mcdc) - else: - # Fission source - sample_phasespace_fission_nuclide(P, nuclide, P_new, mcdc) - - # Assign sensitivity_ID - P_new["sensitivity_ID"] = ID - - # Put the current particle into the secondary bank - adapt.add_active(P_new, prog) - - # ============================================================================== # Particle tracker # ============================================================================== @@ -3890,29 +3524,6 @@ def allocate_pid(P, mcdc): P["track_pid"] = adapt.global_add(mcdc["particle_track_particle_ID"], 0, 1) -# ============================================================================== -# Derivative Source Method (DSM) -# ============================================================================== - - -@njit -def get_DSM_ID(ID1, ID2, Np): - # First-order sensitivity - if ID1 == 0: - return ID2 - - # Self second-order - if ID1 == ID2: - return Np + ID1 - - # Cross second-order - ID1 -= 1 - ID2 -= 1 - return int( - 2 * Np + (Np * (Np - 1) / 2) - (Np - ID1) * ((Np - ID1) - 1) / 2 + ID2 - ID1 - ) - - # ============================================================================= # Continuous Energy Physics # ============================================================================= diff --git a/mcdc/loop.py b/mcdc/loop.py index a614bc19..a3d9b04b 100644 --- a/mcdc/loop.py +++ b/mcdc/loop.py @@ -498,15 +498,6 @@ def step_particle(P, prog): elif event == EVENT_FISSION: kernel.fission(P, prog) - # Sensitivity quantification for nuclide? - material = mcdc["materials"][P["material_ID"]] - if material["sensitivity"] and ( - P["sensitivity_ID"] == 0 - or mcdc["technique"]["dsm_order"] == 2 - and P["sensitivity_ID"] <= mcdc["setting"]["N_sensitivity"] - ): - kernel.sensitivity_material(P, prog) - # Surface crossing if event & EVENT_SURFACE: kernel.surface_crossing(P, prog) @@ -565,7 +556,6 @@ def generate_precursor_particle(DNP, particle_idx, seed_work, prog): P_new["rng_seed"] = part_seed P_new["alive"] = True P_new["w"] = 1.0 - P_new["sensitivity_ID"] = 0 # Set position P_new["x"] = DNP["x"] diff --git a/mcdc/main.py b/mcdc/main.py index d47f7e91..06b9a0c8 100644 --- a/mcdc/main.py +++ b/mcdc/main.py @@ -428,7 +428,7 @@ def prepare(): N_nuclide = len(input_deck.nuclides) for i in range(N_nuclide): # General data - for name in ["ID", "fissionable", "sensitivity", "sensitivity_ID", "dsm_Np"]: + for name in ["ID", "fissionable"]: copy_field(mcdc["nuclides"][i], input_deck.nuclides[i], name) # MG data @@ -765,12 +765,6 @@ def prepare(): mcdc["technique"]["iqmc"]["score"][name] = value # minimum particle weight iqmc["w_min"] = 1e-13 - # ========================================================================= - # Derivative Source Method - # ========================================================================= - - # Threshold - mcdc["technique"]["dsm_order"] = input_deck.technique["dsm_order"] # ========================================================================= # Variance Deconvolution - UQ diff --git a/mcdc/type_.py b/mcdc/type_.py index 43b0e333..72592ca3 100644 --- a/mcdc/type_.py +++ b/mcdc/type_.py @@ -184,7 +184,6 @@ def make_type_particle(input_deck): ("surface_ID", int64), ("translation", float64, (3,)), ("event", int64), - ("sensitivity_ID", int64), ("rng_seed", uint64), ] @@ -223,7 +222,6 @@ def make_type_particle_record(input_deck): ("g", uint64), ("E", float64), ("w", float64), - ("sensitivity_ID", int64), ("rng_seed", uint64), ] @@ -348,9 +346,6 @@ def make_type_nuclide(input_deck): struct = [ ("ID", int64), ("fissionable", bool_), - ("sensitivity", bool_), - ("sensitivity_ID", int64), - ("dsm_Np", float64), ("uq", bool_), ] @@ -446,7 +441,6 @@ def make_type_material(input_deck): struct = [ ("ID", int64), ("N_nuclide", int64), - ("sensitivity", bool_), ("nuclide_IDs", int64, (Nmax_nuclide,)), ("nuclide_densities", float64, (Nmax_nuclide,)), ("uq", bool_), @@ -506,9 +500,6 @@ def make_type_surface(input_deck): ("nx", float64), ("ny", float64), ("nz", float64), - ("sensitivity", bool_), - ("sensitivity_ID", int64), - ("dsm_Np", float64), ] ) @@ -659,14 +650,6 @@ def make_type_source(input_deck): def make_type_tally(input_deck): global tally - # Number of sensitivitys parameters - N_sensitivity = input_deck.setting["N_sensitivity"] - - # Number of tally scores - Ns = 1 + N_sensitivity - if input_deck.technique["dsm_order"] == 2: - Ns = 1 + 2 * N_sensitivity + int(0.5 * N_sensitivity * (N_sensitivity - 1)) - # Get card card = input_deck.tally @@ -688,13 +671,13 @@ def make_type_score(shape): # Scores and shapes scores_shapes = [ - ["flux", (Ns, Ng, Nt, Nx, Ny, Nz, Nmu, N_azi)], - ["density", (Ns, Ng, Nt, Nx, Ny, Nz, Nmu, N_azi)], - ["fission", (Ns, Ng, Nt, Nx, Ny, Nz, Nmu, N_azi)], - ["total", (Ns, Ng, Nt, Nx, Ny, Nz, Nmu, N_azi)], - ["current", (Ns, Ng, Nt, Nx, Ny, Nz, 3)], - ["eddington", (Ns, Ng, Nt, Nx, Ny, Nz, 6)], - ["exit", (Ns, Ng, Nt, 2, Ny, Nz, Nmu, N_azi)], + ["flux", (Ng, Nt, Nx, Ny, Nz, Nmu, N_azi)], + ["density", (Ng, Nt, Nx, Ny, Nz, Nmu, N_azi)], + ["fission", (Ng, Nt, Nx, Ny, Nz, Nmu, N_azi)], + ["total", (Ng, Nt, Nx, Ny, Nz, Nmu, N_azi)], + ["current", (Ng, Nt, Nx, Ny, Nz, 3)], + ["eddington", (Ng, Nt, Nx, Ny, Nz, 6)], + ["exit", (Ng, Nt, 2, Ny, Nz, Nmu, N_azi)], ] # Add score flags to structure @@ -760,8 +743,6 @@ def make_type_setting(deck): ("IC_file", bool_), ("IC_file_name", str_), ("N_precursor", uint64), - # TODO: Move to technique - ("N_sensitivity", uint64), ] # Finalize setting type @@ -990,14 +971,6 @@ def make_type_technique(input_deck): ("IC_fission", float64), ] - # ========================================================================= - # Derivative Source Method - # ========================================================================= - - struct += [ - ("dsm_order", int64), - ] - # ========================================================================= # Variance Deconvolution # ========================================================================= @@ -1022,9 +995,6 @@ def make_type_uq_score(shape): # Tally estimator flags struct = [] - # Number of tally scores - Ns = 1 + input_deck.setting["N_sensitivity"] - # Tally card tally_card = input_deck.tally @@ -1033,13 +1003,13 @@ def make_type_uq_score(shape): # Scores and shapes scores_shapes = [ - ["flux", (Ns, Ng, Nt, Nx, Ny, Nz, Nmu, N_azi)], - ["density", (Ns, Ng, Nt, Nx, Ny, Nz, Nmu, N_azi)], - ["fission", (Ns, Ng, Nt, Nx, Ny, Nz, Nmu, N_azi)], - ["total", (Ns, Ng, Nt, Nx, Ny, Nz, Nmu, N_azi)], - ["current", (Ns, Ng, Nt, Nx, Ny, Nz, 3)], - ["eddington", (Ns, Ng, Nt, Nx, Ny, Nz, 6)], - ["exit", (Ns, Ng, Nt, 2, Ny, Nz, Nmu, N_azi)], + ["flux", (Ng, Nt, Nx, Ny, Nz, Nmu, N_azi)], + ["density", (Ng, Nt, Nx, Ny, Nz, Nmu, N_azi)], + ["fission", (Ng, Nt, Nx, Ny, Nz, Nmu, N_azi)], + ["total", (Ng, Nt, Nx, Ny, Nz, Nmu, N_azi)], + ["current", (Ng, Nt, Nx, Ny, Nz, 3)], + ["eddington", (Ng, Nt, Nx, Ny, Nz, 6)], + ["exit", (Ng, Nt, 2, Ny, Nz, Nmu, N_azi)], ] # Add score flags to structure diff --git a/test/regression/dsm_azurv1/answer.h5 b/test/regression/dsm_azurv1/answer.h5 deleted file mode 100644 index 0c1afb23..00000000 Binary files a/test/regression/dsm_azurv1/answer.h5 and /dev/null differ diff --git a/test/regression/dsm_azurv1/input.py b/test/regression/dsm_azurv1/input.py deleted file mode 100644 index 2b358100..00000000 --- a/test/regression/dsm_azurv1/input.py +++ /dev/null @@ -1,37 +0,0 @@ -import numpy as np -import h5py - -import mcdc - - -# ========================================================================= -# Set model and run -# ========================================================================= - -n1 = mcdc.nuclide(capture=np.array([0.5])) -n2 = mcdc.nuclide( - capture=np.array([0.1]), - fission=np.array([0.4]), - nu_p=np.array([2.5]), - sensitivity=True, -) - -m = mcdc.material(nuclides=[(n1, 1.0), (n2, 1.0)]) - -s1 = mcdc.surface("plane-x", x=-1e10, bc="reflective") -s2 = mcdc.surface("plane-x", x=1e10, bc="reflective") - -mcdc.cell(+s1 & -s2, m) - -mcdc.source(point=[0.0, 0.0, 0.0], isotropic=True) - -scores = ["flux"] -mcdc.tally( - scores=scores, - x=np.linspace(-20.0, 20.0, 202), - t=np.linspace(0.0, 20.0, 21), -) - -mcdc.setting(N_particle=30) - -mcdc.run() diff --git a/test/regression/dsm_lattice/answer.h5 b/test/regression/dsm_lattice/answer.h5 deleted file mode 100644 index ea8ee263..00000000 Binary files a/test/regression/dsm_lattice/answer.h5 and /dev/null differ diff --git a/test/regression/dsm_lattice/c5g7.h5 b/test/regression/dsm_lattice/c5g7.h5 deleted file mode 100644 index 218b300f..00000000 Binary files a/test/regression/dsm_lattice/c5g7.h5 and /dev/null differ diff --git a/test/regression/dsm_lattice/input.py b/test/regression/dsm_lattice/input.py deleted file mode 100644 index 7fc70b08..00000000 --- a/test/regression/dsm_lattice/input.py +++ /dev/null @@ -1,49 +0,0 @@ -import numpy as np -import h5py - -import mcdc - - -# ========================================================================= -# Set model and run -# ========================================================================= - -lib = h5py.File("c5g7.h5", "r") - - -def set_mat(mat): - return mcdc.material( - capture=mat["capture"][:], - scatter=mat["scatter"][:], - fission=mat["fission"][:], - nu_p=mat["nu_p"][:], - nu_d=mat["nu_d"][:], - chi_p=mat["chi_p"][:], - chi_d=mat["chi_d"][:], - speed=mat["speed"], - decay=mat["decay"], - sensitivity=True, - ) - - -mat_uo2 = set_mat(lib["uo2"]) # Fuel: UO2 -mat_mod = set_mat(lib["mod"]) # Moderator -mat_cr = set_mat(lib["cr"]) # Control rod - -s1 = mcdc.surface("plane-x", x=0.0, bc="reflective") -s2 = mcdc.surface("plane-x", x=0.5, sensitivity=True) -s3 = mcdc.surface("plane-x", x=1.5, sensitivity=True) -s4 = mcdc.surface("plane-x", x=2.0, bc="reflective") - -mcdc.cell(+s1 & -s2, mat_uo2) -mcdc.cell(+s2 & -s3, mat_mod) -mcdc.cell(+s3 & -s4, mat_cr) - -mcdc.source(point=[1.0, 0.0, 0.0], energy=[1, 0, 0, 0, 0, 0, 0], isotropic=True) - -scores = ["flux"] -mcdc.tally(scores=scores, x=np.linspace(0.0, 2.0, 11), g="all") - -mcdc.setting(N_particle=3, active_bank_buff=1000) - -mcdc.run() diff --git a/test/unit/test_input_.py b/test/unit/test_input_.py index 59158962..b7c8660e 100644 --- a/test/unit/test_input_.py +++ b/test/unit/test_input_.py @@ -50,7 +50,6 @@ def test_nuclide_basic(): ), speed=np.array([1.0, 2.0, 3.0, 4.0]), decay=np.array([1.0, 2.0, 3.0]), - sensitivity=True, ) assert n1.tag == "Nuclide" @@ -80,8 +79,6 @@ def test_nuclide_basic(): assert (n1.chi_s == np.ones([4, 4]) * 0.25).all() assert (n1.chi_p == np.ones([4, 4]) * 0.25).all() assert (n1.chi_d == np.ones([3, 4]) * 0.25).all() - assert n1.sensitivity == True - assert n1.sensitivity_ID == 1 def test_nuclide_default(): @@ -109,8 +106,6 @@ def test_nuclide_default(): assert (n1.chi_s == np.zeros((5, 5))).all() assert (n1.chi_p == np.zeros((5, 5))).all() assert (n1.chi_d == np.zeros((0, 5))).all() - assert n1.sensitivity == False - assert n1.sensitivity_ID == 0 assert (n2.capture == np.zeros(5)).all() assert (n2.scatter == np.ones(5) * 5.0).all() @@ -122,8 +117,6 @@ def test_nuclide_default(): assert (n2.chi_s == np.ones((5, 5)) * 0.2).all() assert (n2.chi_p == np.zeros((5, 5))).all() assert (n2.chi_d == np.zeros((0, 5))).all() - assert n2.sensitivity == False - assert n2.sensitivity_ID == 0 assert (n3.capture == np.zeros(5)).all() assert (n3.scatter == np.zeros(5)).all() @@ -135,45 +128,6 @@ def test_nuclide_default(): assert (n3.chi_s == np.zeros((5, 5))).all() assert (n3.chi_p == np.ones((5, 5)) * 0.2).all() assert (n3.chi_d == np.zeros((0, 5))).all() - assert n3.sensitivity == False - assert n3.sensitivity_ID == 0 - - -def test_nuclide_IDs_sensitivity(): - # Start fresh - mcdc.reset() - - # Create nuclides with various sensitivity tags - n1 = mcdc.nuclide(capture=np.array([1.0]), sensitivity=True) - n2 = mcdc.nuclide(capture=np.array([1.0])) - n3 = mcdc.nuclide(capture=np.array([1.0]), sensitivity=True) - n4 = mcdc.nuclide(capture=np.array([1.0]), sensitivity=True) - n5 = mcdc.nuclide(capture=np.array([1.0])) - n6 = mcdc.nuclide(capture=np.array([1.0])) - n7 = mcdc.nuclide(capture=np.array([1.0]), sensitivity=True) - - # Checks - assert n1.ID == 0 - assert n2.ID == 1 - assert n3.ID == 2 - assert n4.ID == 3 - assert n5.ID == 4 - assert n6.ID == 5 - assert n7.ID == 6 - assert n1.sensitivity - assert not n2.sensitivity - assert n3.sensitivity - assert n4.sensitivity - assert not n5.sensitivity - assert not n6.sensitivity - assert n7.sensitivity - assert n1.sensitivity_ID == 1 - assert n2.sensitivity_ID == 0 - assert n3.sensitivity_ID == 2 - assert n4.sensitivity_ID == 3 - assert n5.sensitivity_ID == 0 - assert n6.sensitivity_ID == 0 - assert n7.sensitivity_ID == 4 # ====================================================================================== @@ -224,7 +178,6 @@ def test_material_single(): ), speed=np.array([1.0, 2.0, 3.0, 4.0]), decay=np.array([1.0, 2.0, 3.0]), - sensitivity=True, ) # Checks @@ -256,7 +209,6 @@ def test_material_single(): assert (m1.nu_f == np.array([4.0, 8.0, 12.0, 16.0])).all() assert (m1.chi_s == np.ones([4, 4]) * 0.25).all() assert (m1.chi_p == np.ones([4, 4]) * 0.25).all() - assert m1.sensitivity == True # Check if the nuclide was registered n2 = mcdc.nuclide(capture=np.ones(5)) @@ -269,7 +221,7 @@ def test_material_multi(): # Create a multi-nuclide material n1 = mcdc.nuclide(capture=np.ones(5), speed=np.ones(5) * 1) - n2 = mcdc.nuclide(scatter=np.ones((5, 5)), speed=np.ones(5) * 2, sensitivity=True) + n2 = mcdc.nuclide(scatter=np.ones((5, 5)), speed=np.ones(5) * 2) n3 = mcdc.nuclide( fission=np.ones(5), nu_p=np.ones(5), chi_p=np.ones((5, 5)), speed=np.ones(5) * 3 ) @@ -285,37 +237,6 @@ def test_material_multi(): assert m1.J == 0 assert (m1.total == np.ones(5) * 14).all() assert (m1.speed == np.ones(5) * 30 / 14).all() - assert m1.sensitivity == True - - -def test_material_IDs_sensitivity(): - # Start fresh - mcdc.reset() - - # Create materials with various sensitivity tags - m1 = mcdc.material(capture=np.array([1.0]), sensitivity=True) - m2 = mcdc.material(capture=np.array([1.0])) - m3 = mcdc.material(capture=np.array([1.0]), sensitivity=True) - m4 = mcdc.material(capture=np.array([1.0]), sensitivity=True) - m5 = mcdc.material(capture=np.array([1.0])) - m6 = mcdc.material(capture=np.array([1.0])) - m7 = mcdc.material(capture=np.array([1.0]), sensitivity=True) - - # Checks - assert m1.ID == 0 - assert m2.ID == 1 - assert m3.ID == 2 - assert m4.ID == 3 - assert m5.ID == 4 - assert m6.ID == 5 - assert m7.ID == 6 - assert m1.sensitivity - assert not m2.sensitivity - assert m3.sensitivity - assert m4.sensitivity - assert not m5.sensitivity - assert not m6.sensitivity - assert m7.sensitivity # ======================================================================================