From a258658f130921a235754f6eda44462eef0400ac Mon Sep 17 00:00:00 2001 From: "Joanna Piper Morgan (jonsey)" Date: Thu, 24 Apr 2025 16:02:43 -0700 Subject: [PATCH 01/17] adding delta tracking function for input decks --- mcdc/__init__.py | 1 + 1 file changed, 1 insertion(+) diff --git a/mcdc/__init__.py b/mcdc/__init__.py index 6851a6f4..ec7a9eaf 100644 --- a/mcdc/__init__.py +++ b/mcdc/__init__.py @@ -20,6 +20,7 @@ weight_roulette, IC_generator, uq, + delta_tracking, reset, domain_decomposition, make_particle_bank, From e9f6456fb929c8f8e61a79abbb7a89ce54984990 Mon Sep 17 00:00:00 2001 From: "Joanna Piper Morgan (jonsey)" Date: Thu, 24 Apr 2025 16:04:14 -0700 Subject: [PATCH 02/17] adding delta tracking types including the majorant into technique, also keeping track of how many boundary surfaces there are --- mcdc/type_.py | 42 ++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 42 insertions(+) diff --git a/mcdc/type_.py b/mcdc/type_.py index d005d057..8a8b4ded 100644 --- a/mcdc/type_.py +++ b/mcdc/type_.py @@ -221,6 +221,7 @@ def make_type_particle(input_deck): ("fresh", bool_), ("event", int64), ("rng_seed", uint64), + ("delta_tracking", bool_), ] # Get modes @@ -1100,6 +1101,10 @@ def make_type_technique(input_deck): ("branchless_collision", bool_), ("domain_decomposition", bool_), ("uq", bool_), + ("delta_tracking", bool_), + ("collision_estimator", bool_), + ("dt_cutoff_energy", float64), # cut off energy for delta tracking + ("dt_material_exclusion", int64), # exculsion materials for delta tracking ] # ========================================================================= @@ -1287,10 +1292,38 @@ def make_type_technique(input_deck): struct += [("uq_", uq)] + # ========================================================================= + # Delta Tracking + # ========================================================================= + + if mode_MG: + struct += [("majorant_energy", float64, (G,))] + struct += [("micro_majorant_xsec", float64, (G,))] + else: + N_majorant = compute_majorant_size(input_deck) + struct += [("majorant_energy", float64, (N_majorant,))] + struct += [("micro_majorant_xsec", float64, (N_majorant,))] + + struct += [("N_majorant", int64)] + # Finalize technique type technique = into_dtype(struct) +def compute_majorant_size(input_deck): + """just computes size of majorant! + this will be redone in prepare()""" + + energy_grid = [] + + dir_name = os.getenv("MCDC_XSLIB") + for nuc in input_deck.nuclides: + with h5py.File(dir_name + "/" + nuc.name + ".h5", "r") as f: + energy_grid = np.append(energy_grid, f["E_xs"][:]) + + return np.size(np.unique(energy_grid)) + 1 + + # UQ def make_type_uq(input_deck): global uq, uq_nuc, uq_mat @@ -1449,6 +1482,14 @@ def make_type_global(input_deck): N_cell_tally = len(input_deck.cell_tallies) N_cs_tally = len(input_deck.cs_tallies) + # find the number of reflecting surfaces + N_boundary = 0 + for i in range(N_surface): + if input_deck.surfaces[i].boundary_type == "reflective": + N_boundary += 1 + elif input_deck.surfaces[i].boundary_type == "vacuum": + N_boundary += 1 + # Cell data sizes N_cell_surface = sum([len(x.surface_IDs) for x in input_deck.cells]) N_cell_region = sum([len(x._region_RPN) for x in input_deck.cells]) @@ -1515,6 +1556,7 @@ def make_type_global(input_deck): ("nuclides", nuclide, (N_nuclide,)), ("materials", material, (N_material,)), ("surfaces", surface, (N_surface,)), + ("boundary_surface_IDs", int64, (N_boundary,)), # delta tracking # Cells ("cells", cell, (N_cell,)), ("cells_data_surface", int64, (N_cell_surface,)), From b09e705a43645fcdf1e36ad019fe9d44485497e9 Mon Sep 17 00:00:00 2001 From: "Joanna Piper Morgan (jonsey)" Date: Thu, 24 Apr 2025 16:05:55 -0700 Subject: [PATCH 03/17] command line options for delta tracking --- mcdc/config.py | 7 +++++++ 1 file changed, 7 insertions(+) diff --git a/mcdc/config.py b/mcdc/config.py index 3b80f3d0..4877efad 100644 --- a/mcdc/config.py +++ b/mcdc/config.py @@ -67,6 +67,12 @@ parser.add_argument("--caching", action="store_true") parser.add_argument("--no_caching", dest="caching", action="store_false") parser.add_argument("--runtime_output", default=False, action="store_true") +parser.add_argument( + "--delta_tracking", + default=False, + action="store_true", + help="Set delta tracking to on", +) parser.set_defaults(caching=False) args, unargs = parser.parse_known_args() @@ -77,6 +83,7 @@ clear_cache = args.clear_cache from mpi4py import MPI + import shutil src_path = os.path.dirname(os.path.abspath(__file__)) From e4fae1286930e0282f7f975c223b9f7d6db94f63 Mon Sep 17 00:00:00 2001 From: "Joanna Piper Morgan (jonsey)" Date: Thu, 24 Apr 2025 16:06:32 -0700 Subject: [PATCH 04/17] event_phantom_collision and majorant declerations --- mcdc/constant.py | 2 ++ 1 file changed, 2 insertions(+) diff --git a/mcdc/constant.py b/mcdc/constant.py index 1667f65c..ec96afe4 100644 --- a/mcdc/constant.py +++ b/mcdc/constant.py @@ -83,6 +83,7 @@ EVENT_TIME_CENSUS = 1 << 9 EVENT_TIME_BOUNDARY = 1 << 10 EVENT_IQMC_MESH = 1 << 11 +EVENT_PHANTOM_COLLISION = 1 << 12 # Gyration raius type GYRATION_RADIUS_ALL = 0 @@ -148,6 +149,7 @@ XS_NU_FISSION_PROMPT = 5 XS_NU_FISSION_DELAYED = 6 XS_NU_SCATTER = 7 +XS_MAJORANT = 8 NU_FISSION = 0 NU_FISSION_PROMPT = 1 From e3b0f78b0b05a3c22e82586902f43158993b01f6 Mon Sep 17 00:00:00 2001 From: "Joanna Piper Morgan (jonsey)" Date: Thu, 24 Apr 2025 16:07:08 -0700 Subject: [PATCH 05/17] delta tracking input options including plotting the majorand energy cutoff and boolian operators --- mcdc/global_.py | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/mcdc/global_.py b/mcdc/global_.py index 64895cca..f062df31 100644 --- a/mcdc/global_.py +++ b/mcdc/global_.py @@ -153,6 +153,11 @@ def reset(self): "IC_cycle_stretch": 1.0, "branchless_collision": False, "uq": False, + "delta_tracking": False, + "collision_estimator": False, + "dt_cutoff_energy": 0, + "dt_material_exclusion": INF, + "plot_majorant": True, } self.uq_deltas = { From 4011af24dc117f0c4b3b1ed90e8c9003a41649d9 Mon Sep 17 00:00:00 2001 From: "Joanna Piper Morgan (jonsey)" Date: Thu, 24 Apr 2025 16:07:28 -0700 Subject: [PATCH 06/17] delta tracking input function --- mcdc/input_.py | 44 ++++++++++++++++++++++++++++++++++++++++---- 1 file changed, 40 insertions(+), 4 deletions(-) diff --git a/mcdc/input_.py b/mcdc/input_.py index fa29af35..080ef977 100644 --- a/mcdc/input_.py +++ b/mcdc/input_.py @@ -227,6 +227,7 @@ def material( chi_d=None, speed=None, decay=None, + name=None, ): """ Create a material @@ -258,6 +259,8 @@ def material( Energy group speed [cm/s]. decay : numpy.ndarray (1D), optional Precursor group decay constant [/s]. + name : str, optional + Name of the material (e.g. moderator) Returns ------- @@ -324,7 +327,7 @@ def material( print_error( "Continuous energy data directory not configured \n " "see https://cement-psaapgithubio.readthedocs.io/en/latest" - "/install.html#configuring-continuous-energy-library \n" + "./install.html#configuring-continuous-energy-library \n" ) # Fissionable flag @@ -446,8 +449,7 @@ def material( def surface(type_, bc="interface", **kw): - """ - Create a surface to define the region of a cell. + """Create a surface to define the region of a cell. Parameters ---------- @@ -1245,6 +1247,40 @@ def branchless_collision(): card["weighted_emission"] = False +def delta_tracking( + collision_estimator=False, + cutoff_energy=0, + material_exclusion=INF, + plot_majorant=True, +): + """ + Activate delta tracking and set hybrid parameters + + Parameters + ---------- + collision_estimator : bool + Use a collision estimator for flux tallies + cutoff_energy : float + A cut off energy, underwich do surface tracking (**continuous energy only**) (look at the majorant plot and make an informed decision!) + material_exclusion : int + Material ID of a single region not to do delta tracking in (TODO: Make multiple regions) (**continuous energy only**) + plot_majorant : bool + Saves a .pdf plot of the majorant and nuclides when produced + + Notes + ----- + TODO: turn off other incompatible tracking techniques if any (k-eigenvalue, fission flux tallies, ) + Doesn't work with k-eigenvalue simulations + + """ + card = global_.input_deck.technique + card["delta_tracking"] = True + card["collision_estimator"] = collision_estimator + card["dt_cutoff_energy"] = cutoff_energy + card["dt_material_exclusion"] = material_exclusion + card["plot_majorant"] = plot_majorant + + def time_census(t, tally_frequency=None): """ Set time-census boundaries. @@ -1253,7 +1289,7 @@ def time_census(t, tally_frequency=None): ---------- t : array_like[float] The time-census boundaries. - tally_frecuency : integer, optional + tally_frequency : integer, optional Number of uniform tally time mesh bins in census-based tallying. This overrides manual tally time mesh definitions. From 88669a74e33bff21a91dd0fd128813638c0f02a7 Mon Sep 17 00:00:00 2001 From: "Joanna Piper Morgan (jonsey)" Date: Thu, 24 Apr 2025 16:08:15 -0700 Subject: [PATCH 07/17] hybrid delta tracking, collision estimator, rejection samepling, majorant xsec lookups --- mcdc/kernel.py | 240 +++++++++++++++++++++++++++++++++++++++++++------ 1 file changed, 214 insertions(+), 26 deletions(-) diff --git a/mcdc/kernel.py b/mcdc/kernel.py index 8305716e..06e645ce 100644 --- a/mcdc/kernel.py +++ b/mcdc/kernel.py @@ -1849,6 +1849,79 @@ def mesh_get_energy_index(P_arr, mesh, mode_MG): # ============================================================================= +@njit +def score_mesh_collision_tally(P_arr, tally, data, mcdc): + """Scoring tallies with a collision estimator,""" + P = P_arr[0] + tally_bin = data[TALLY] + material = mcdc["materials"][P["material_ID"]] + mesh = tally["filter"] + stride = tally["stride"] + + # Easily identified tally bin indices + mu, azi = mesh_get_angular_index(P_arr, mesh) + g, outside_energy = mesh_get_energy_index(P_arr, mesh, mcdc["setting"]["mode_MG"]) + + # Get current indices indices + ix, iy, iz, it, outside = mesh_.structured.get_indices(P_arr, mesh) + + if outside or outside_energy: + return + + idx = ( + stride["tally"] + + mu * stride["mu"] + + azi * stride["azi"] + + g * stride["g"] + + it * stride["t"] + + ix * stride["x"] + + iy * stride["y"] + + iz * stride["z"] + ) + + # collision estimator requires all + SigmaT = get_MacroXS(XS_TOTAL, material, P_arr, mcdc) + flux = 1 / SigmaT + + mu = P["ux"] + for i in range(tally["N_score"]): + score_type = tally["scores"][i] + score = 0 + if score_type == SCORE_FLUX: + score = flux + elif score_type == SCORE_DENSITY: + ut = 1.0 / physics.get_speed(P_arr, mcdc) + score = flux * ut + elif score_type == SCORE_TOTAL: + score = flux * SigmaT + elif score_type == SCORE_FISSION: + SigmaF = get_MacroXS(XS_FISSION, material, P_arr, mcdc) + score = flux * SigmaF + if score_type == SCORE_NET_CURRENT: + score = flux * mu + if score_type == SCORE_MU_SQ: + score = flux * mu * mu + # elif score_type == SCORE_TIME_MOMENT_FLUX: + # score = flux * (t - (mesh["t"][it - 1] + mesh["t"][it]) / 2) + # elif score_type == SCORE_SPACE_MOMENT_FLUX: + # score = flux * (x - (mesh["x"][ix + 1] + mesh["x"][ix]) / 2) + # elif score_type == SCORE_TIME_MOMENT_CURRENT: + # score = flux * mu * (t - (mesh["t"][it - 1] + mesh["t"][it]) / 2) + # elif score_type == SCORE_SPACE_MOMENT_CURRENT: + # score = flux * mu * (x - (mesh["x"][ix + 1] + mesh["x"][ix]) / 2) + # elif score_type == SCORE_TIME_MOMENT_MU_SQ: + # score = flux * mu * mu * (t - (mesh["t"][it - 1] + mesh["t"][it]) / 2) + # elif score_type == SCORE_SPACE_MOMENT_MU_SQ: + # score = flux * mu * mu * (x - (mesh["x"][ix + 1] + mesh["x"][ix]) / 2) + adapt.global_add(tally_bin, (TALLY_SCORE, idx + i), round(score)) + + SigmaT = get_MacroXS(XS_TOTAL, material, P_arr, mcdc) + score = 1.0 / SigmaT + adapt.global_add(tally_bin, (TALLY_SCORE, idx + i), round(score)) + + return + + @njit def score_mesh_tally(P_arr, distance, tally, data, mcdc): P = P_arr[0] @@ -2768,11 +2841,33 @@ def move_to_event(P_arr, data, mcdc): # ================================================================================== P = P_arr[0] - # Multigroup preparation - # In MG mode, particle speed is material-dependent. - if mcdc["setting"]["mode_MG"]: - # If material is not identified yet, locate the particle - if P["material_ID"] == -1: + # Delta tracking logic + # TODO: should be it's own function of logic, maybe setting in loop + if mcdc["technique"]["delta_tracking"]: + if mcdc["setting"]["mode_MG"]: + P["delta_tracking"] = True + else: + P["delta_tracking"] = True + + # energy cutoff (above delta tracking, below surface) + if P["E"] < mcdc["technique"]["dt_cutoff_energy"]: + P["delta_tracking"] = False + + # don't delta track if in excluded material + if P["material_ID"] == mcdc["technique"]["dt_material_exclusion"]: + P["delta_tracking"] = False + + # but complete rejection sampleing process before switching + if P["event"] & EVENT_PHANTOM_COLLISION: + P["delta_tracking"] = True + else: + P["delta_tracking"] = False + + # If material is not identified yet, locate the particle + if P["material_ID"] == -1: + # Multigroup preparation + # In MG mode, particle speed is material-dependent. + if mcdc["setting"]["mode_MG"]: if not geometry.locate_particle(P_arr, mcdc): # Particle is lost P["event"] = EVENT_LOST @@ -2786,12 +2881,19 @@ def move_to_event(P_arr, data, mcdc): # - Set particle boundary event (surface or lattice crossing, or lost) # - Return distance to boundary (surface or lattice) - d_boundary = geometry.inspect_geometry(P_arr, mcdc) + if P["delta_tracking"]: + d_boundary = geometry.distance_to_boundary(P_arr, mcdc) + else: + d_boundary = geometry.inspect_geometry(P_arr, mcdc) # Particle is lost? if P["event"] == EVENT_LOST: return + # Particle is dead? + if P["alive"] == False: + return + # ================================================================================== # Get distances to other events # ================================================================================== @@ -2851,40 +2953,65 @@ def move_to_event(P_arr, data, mcdc): P["event"] = EVENT_TIME_BOUNDARY P["surface_ID"] = -1 + # Score tracklength tallies + if not mcdc["technique"]["collision_estimator"]: + if mcdc["cycle_active"]: + # Mesh tallies + for tally in mcdc["mesh_tallies"]: + score_mesh_tally(P_arr, distance, tally, data, mcdc) + + # Cell tallies + cell = mcdc["cells"][P["cell_ID"]] + for i in range(cell["N_tally"]): + ID = cell["tally_IDs"][i] + tally = mcdc["cell_tallies"][ID] + score_cell_tally(P_arr, distance, tally, data, mcdc) + + # CS tallies + for tally in mcdc["cs_tallies"]: + score_cs_tally(P_arr, distance, tally, data, mcdc) + + if mcdc["setting"]["mode_eigenvalue"] and not P["delta_tracking"]: + eigenvalue_tally(P_arr, distance, mcdc) + # ========================================================================= # Move particle # ========================================================================= - # Score tracklength tallies - if mcdc["cycle_active"]: - # Mesh tallies - for tally in mcdc["mesh_tallies"]: - score_mesh_tally(P_arr, distance, tally, data, mcdc) + move_particle(P_arr, distance, mcdc) - # Cell tallies - cell = mcdc["cells"][P["cell_ID"]] - for i in range(cell["N_tally"]): - ID = cell["tally_IDs"][i] - tally = mcdc["cell_tallies"][ID] - score_cell_tally(P_arr, distance, tally, data, mcdc) + # ========================================================================= + # Delta Tracking + # ========================================================================= - # CS tallies - for tally in mcdc["cs_tallies"]: - score_cs_tally(P_arr, distance, tally, data, mcdc) + # Delta tracking rejection sampleing must happen AFTER the particle is moved + if P["delta_tracking"] and P["event"] == EVENT_COLLISION: + # finding the particle + P["cell_ID"] = -1 + P["material_ID"] = -1 + geometry.locate_particle(P_arr, mcdc) - if mcdc["setting"]["mode_eigenvalue"]: - eigenvalue_tally(P_arr, distance, mcdc) + rejection_sample(P_arr, mcdc) - # Move particle - move_particle(P_arr, distance, mcdc) + if mcdc["technique"]["collision_estimator"]: + if P["event"] == EVENT_COLLISION and mcdc["cycle_active"]: + # Mesh tallies + # use a collision estimator for currently undefined tallies + for tally in mcdc["mesh_tallies"]: + if tally is not SCORE_FLUX: + score_mesh_collision_tally(P_arr, tally, data, mcdc) @njit def distance_to_collision(P_arr, mcdc): P = P_arr[0] # Get total cross-section - material = mcdc["materials"][P["material_ID"]] - SigmaT = get_MacroXS(XS_TOTAL, material, P_arr, mcdc) + + if P["delta_tracking"]: + SigmaT = get_MacroMaj(P_arr, mcdc) + else: + material = mcdc["materials"][P["material_ID"]] + SigmaT = get_MacroXS(XS_TOTAL, material, P_arr, mcdc) # Vacuum material? if SigmaT == 0.0: @@ -2893,9 +3020,29 @@ def distance_to_collision(P_arr, mcdc): # Sample collision distance xi = rng(P_arr) distance = -math.log(xi) / SigmaT + return distance +@njit +def rejection_sample(P_arr, mcdc): + # For use in delta tracking + P = P_arr[0] + + material = mcdc["materials"][P["material_ID"]] + SigmaT = get_MacroXS(XS_TOTAL, material, P_arr, mcdc) + Maj = get_MacroMaj(P_arr, mcdc) + + reject_rat = SigmaT / Maj + + xi = rng(P_arr) + + if xi < reject_rat: # real collision + P["event"] = EVENT_COLLISION + else: # phantom collision + P["event"] = EVENT_PHANTOM_COLLISION + + # ============================================================================= # Surface crossing # ============================================================================= @@ -3710,6 +3857,46 @@ def weight_roulette(P_arr, mcdc): P["alive"] = False +# ============================================================================= +# Majorant function --- Delta Tracking +# ============================================================================= + + +@njit +def get_MacroMaj(P_arr, mcdc): + """Hybrid delta tracking""" + + P = P_arr[0] + + if mcdc["setting"]["mode_MG"]: + return mcdc["technique"]["micro_majorant_xsec"][P["g"]] + + elif mcdc["setting"]["mode_CE"]: + # nuclide density over all nuclides + # same as assuming all atoms are of same xsec + # N = 0.0 + # for i in range(material["N_nuclide"]): + # N += material["nuclide_densities"][i] + + data = mcdc["technique"]["micro_majorant_xsec"] + Egrid = mcdc["technique"]["majorant_energy"] + E = P["E"] + + Ne = mcdc["technique"]["N_majorant"] + + MacroXS = get_XS(data, E, Egrid, Ne) + + # import matplotlib.pyplot as plt + # plt.figure() + # plt.loglog(Egrid, data) + # plt.loglog(P["E"], MacroXS, "X") + # plt.show() + + # MacroXS *= N + + return MacroXS + + # ============================================================================= # Continuous Energy Physics # ============================================================================= @@ -3757,6 +3944,7 @@ def get_MacroXS(type_, material, P_arr, mcdc): # Sum over all nuclides for i in range(material["N_nuclide"]): + ID_nuclide = material["nuclide_IDs"][i] nuclide = mcdc["nuclides"][ID_nuclide] From 74b857710cacbed6dc36a616c512ca0ea833b511 Mon Sep 17 00:00:00 2001 From: "Joanna Piper Morgan (jonsey)" Date: Thu, 24 Apr 2025 16:10:02 -0700 Subject: [PATCH 08/17] adding a phantom virtual delta tracking operation --- mcdc/loop.py | 7 +++++-- 1 file changed, 5 insertions(+), 2 deletions(-) diff --git a/mcdc/loop.py b/mcdc/loop.py index 6a603026..36acafb6 100644 --- a/mcdc/loop.py +++ b/mcdc/loop.py @@ -434,7 +434,6 @@ def loop_source(seed, data, mcdc): for idx_work in range(work_size): generate_source_particle(work_start, idx_work, seed, mcdc) - # Run the source particle and its secondaries exhaust_active_bank(data, mcdc) source_closeout(mcdc, idx_work, N_prog, data) @@ -566,7 +565,6 @@ def gpu_loop_source(seed, data, mcdc): def loop_particle(P_arr, data, prog): P = P_arr[0] mcdc = adapt.mcdc_global(prog) - while P["alive"]: step_particle(P_arr, data, prog) @@ -581,6 +579,7 @@ def step_particle(P_arr, data, prog): # Execute events if P["event"] == EVENT_LOST: + P["alive"] = False return # Collision @@ -608,6 +607,10 @@ def step_particle(P_arr, data, prog): elif P["event"] & EVENT_FISSION: kernel.fission(P_arr, prog) + if P["event"] & EVENT_PHANTOM_COLLISION: + # delta tracking + P["alive"] = True + # Surface and domain crossing if P["event"] & EVENT_SURFACE_CROSSING: kernel.surface_crossing(P_arr, data, prog) From d67c8fc38159188f43d55eb01b2f7f8c2b48f1ce Mon Sep 17 00:00:00 2001 From: "Joanna Piper Morgan (jonsey)" Date: Thu, 24 Apr 2025 16:10:55 -0700 Subject: [PATCH 09/17] macro and micro scopic continious energy and multigroup majorant computations, visualization improvemtns including a movie maker and legend --- mcdc/main.py | 252 ++++++++++++++++++++++++++++++++++++++++++++++++++- 1 file changed, 250 insertions(+), 2 deletions(-) diff --git a/mcdc/main.py b/mcdc/main.py index a0c96274..9be6e2d1 100644 --- a/mcdc/main.py +++ b/mcdc/main.py @@ -53,7 +53,8 @@ def run(): input_deck.setting["output_name"] = config.args.output if config.args.progress_bar is not None: input_deck.setting["progress_bar"] = config.args.progress_bar - + if config.args.delta_tracking is not None: + input_deck.technique["delta_tracking"] = config.args.delta_tracking # Start timer total_start = MPI.Wtime() @@ -70,7 +71,7 @@ def run(): mcdc["runtime_preparation"] = MPI.Wtime() - preparation_start # Print banner, hardware configuration, and header - print_banner(mcdc) + print_banner(mcdc, config.args) print_msg(" Now running TNT...") if mcdc["setting"]["mode_eigenvalue"]: @@ -644,7 +645,9 @@ def prepare(): # ========================================================================= N_surface = len(input_deck.surfaces) + j = 0 # indexing reflective surfaces for i in range(N_surface): + surface = mcdc["surfaces"][i] surface_input = input_deck.surfaces[i] @@ -694,8 +697,12 @@ def prepare(): mcdc["surfaces"][i]["BC"] = BC_NONE elif input_deck.surfaces[i].boundary_type == "vacuum": mcdc["surfaces"][i]["BC"] = BC_VACUUM + # mcdc["boundary_surface_IDs"][j] = i + # j += 1 # increasing the reflecting surface indexer elif input_deck.surfaces[i].boundary_type == "reflective": mcdc["surfaces"][i]["BC"] = BC_REFLECTIVE + mcdc["boundary_surface_IDs"][j] = i + j += 1 # increasing the reflecting surface indexer # Variables with possible different sizes for name in ["tally_IDs"]: @@ -1227,6 +1234,10 @@ def prepare(): "IC_generator", "branchless_collision", "uq", + "delta_tracking", + "collision_estimator", + "dt_cutoff_energy", + "dt_material_exclusion", ]: copy_field(mcdc["technique"], input_deck.technique, name) @@ -1449,6 +1460,194 @@ def prepare(): if flags["nu_p"] or flags["nu_d"]: flags["nu_f"] = True + # ========================================================================= + # Delta tracking (majorant copmutation) + # ========================================================================= + # TODO: Make seperate function that can be called without run + + micro_maj = False + + if mcdc["technique"]["delta_tracking"]: + + if input_deck.technique["plot_majorant"]: + import matplotlib.pyplot as plt + + plt.figure(figsize=(11, 8.5)) + + if mode_CE: + + # TODO: Use numpy functions (currently can't as data is disorginized and requires extrapolation) + try: + import scipy + except: + print_error( + "SciPy not found!\n" + "Contininous energy delta tracking requires\n\ + interpolation functions from Scipy to construct the marjoant,\n \ + try `pip install scipy`" + ) + + # unify the energy grids from all nuclides + majorant_energy_grid = np.array([]) + for n in range(N_nuclide): + nuclide = mcdc["nuclides"][n] + majorant_energy_grid = np.append(majorant_energy_grid, nuclide["E_xs"]) + + # sort energy grid and eliminate duplicate points + majorant_energy_grid = np.unique(majorant_energy_grid) + majorant_xsec = np.zeros_like(majorant_energy_grid) + + # macro maj + if micro_maj: + + # find the majorant at every point on the energy grid + for n in range(N_nuclide): + nuclide = mcdc["nuclides"][n] + + # builds an interploation funciton for a given nuclide + total_xsec_unified = scipy.interpolate.interp1d( + nuclide["E_xs"], nuclide["ce_total"], bounds_error=False + ) + + # evaluates nuclide interploation function at all unified energy grid points + total_xsec_unified = total_xsec_unified(majorant_energy_grid) + + # compares old majorant xsec and the currently evaluated unified xsec and picks the larger xsecs + majorant_xsec = np.max((majorant_xsec, total_xsec_unified), axis=0) + + import matplotlib.pyplot as plt + + if input_deck.technique["plot_majorant"]: + plt.plot( + majorant_energy_grid[::5], + majorant_xsec[::5], + "k2", + label="majorant", + linewidth=0.75, + ) + for n in range(N_nuclide): + nuclide = mcdc["nuclides"][n] + plt.plot( + nuclide["E_xs"], + nuclide["ce_total"], + label=input_deck.nuclides[n].name, + linewidth=0.75, + ) + plt.xscale("log") + plt.yscale("log") + plt.xlabel("E [ev]") + plt.ylabel("σ") + plt.legend() + plt.savefig("micro_majorant.pdf") + + else: + import matplotlib.pyplot as plt + + plt.figure(figsize=(11, 8.5)) + + for m in range(N_material): + + material = mcdc["materials"][m] + + material_energy_grid = np.array([]) + + # copmute a unified energy grid across all nuclides of a given material + for n in range(material["N_nuclide"]): + nuclide = mcdc["nuclides"][n] + material_energy_grid = np.append( + material_energy_grid, nuclide["E_xs"] + ) + material_energy_grid = np.unique(material_energy_grid) + MacroXS = np.zeros_like(material_energy_grid) + + # compute the macroscopic total cross section of a material on its unified energy grid + for n in range(material["N_nuclide"]): + ID_nuclide = material["nuclide_IDs"][n] + nuclide = mcdc["nuclides"][ID_nuclide] + + # Get nuclide density + N = material["nuclide_densities"][n] + + # putting the microscopic cross-sections on the unifed material energy grid + total_micro_xsec_unified = scipy.interpolate.interp1d( + nuclide["E_xs"], nuclide["ce_total"], bounds_error=False + ) + total_micro_xsec_unified = total_micro_xsec_unified( + material_energy_grid + ) + + # Accumulate + MacroXS += N * total_micro_xsec_unified + + if input_deck.technique["plot_majorant"]: + plt.plot( + material_energy_grid, + MacroXS, + linewidth=0.5, + label="Mat {}".format(m), + ) # label=input_deck.materials[m].name, + + # puting the total macroscopic cross sections on on the majorant energy grid + total_xsec_unified = scipy.interpolate.interp1d( + material_energy_grid, MacroXS, bounds_error=False + ) + total_xsec_unified = total_xsec_unified(majorant_energy_grid) + + # compares old majorant xsec and the currently evaluated unified xsec and picks the larger xsecs + majorant_xsec = np.max((majorant_xsec, total_xsec_unified), axis=0) + + if input_deck.technique["plot_majorant"]: + plt.plot( + majorant_energy_grid[::5], + majorant_xsec[::5], + "k2", + label="majorant", + linewidth=0.5, + ) + plt.xscale("log") + plt.yscale("log") + plt.xlabel("E [kev]") + plt.ylabel("σ") + plt.legend() + plt.savefig("macro_majorant.pdf") + + elif mode_MG: + N_groups = mcdc["nuclides"][0]["G"] + majorant_energy_grid = np.zeros(N_groups) + majorant_xsec = np.zeros(N_groups) + + majorant_energy_grid = mcdc["nuclides"][0]["speed"][:] + + for i in range(N_groups): + # print("g {} ".format(i), mcdc["nuclides"][:]["total"][:, i]) + majorant_xsec[i] = np.max(mcdc["nuclides"][:]["total"][:, i]) + + # print("here!") + import matplotlib.pyplot as plt + + if input_deck.technique["plot_majorant"]: + plt.figure(figsize=(11, 8.5)) + plt.plot(majorant_energy_grid, majorant_xsec, "k2", label="majorant") + for n in range(N_nuclide): + nuclide = mcdc["nuclides"][n] + plt.plot( + nuclide["speed"][:], + nuclide["total"], + "-*", + label=input_deck.nuclides[n].name, + ) + plt.xscale("log") + plt.yscale("log") + plt.xlabel("E [kev]") + plt.ylabel("σ") + plt.legend() + plt.savefig("majorant.pdf") + + # This might need to use copy_feild but not sure why + mcdc["technique"]["micro_majorant_xsec"] = majorant_xsec[:] + mcdc["technique"]["majorant_energy"] = majorant_energy_grid[:] + mcdc["technique"]["N_majorant"] = np.size(majorant_xsec) + # ========================================================================= # MPI # ========================================================================= @@ -2410,7 +2609,10 @@ def visualize( z=0.0, pixels=(100, 100), colors=None, + labels=None, + legend=False, time=np.array([0.0]), + movie=False, save_as=None, ): """ @@ -2432,9 +2634,30 @@ def visualize( Number of respective pixels in the two axes in vis_plane colors : array_like List of pairs of material and its color + labels : array_like + List of pairs of matieral and its label on plot, sets legend to True + legend : bool + Display legend or not, if labels not supplied uses materialIDs + save_as : str + If supplied outputs saved as a time stampped .png instead of displayed + movie : bool + Save output as a .gif, save_as must be supplied """ # TODO: add input error checkers + if movie: + try: + assert save_as is not None + except: + print_error("To save a move supply name of file") + + try: + import imageio + except: + print_error( + "Imageio required to save visualization results as a git\ntry: pip install imageio" + ) + _, mcdc_container = prepare() mcdc = mcdc_container[0] @@ -2450,6 +2673,20 @@ def visualize( colors[i] = plt.cm.Set1(i)[:-1] WHITE = mpl_colors.to_rgb("white") + # legend label assigment (by material ID) + import matplotlib.patches as mpatches + + legend_hands = [0] * len(mcdc["materials"]) + if labels is not None: + legend = True + for item in labels.items(): + legend_hands[item[0].ID] = mpatches.Patch( + color=colors[item[0].ID], label=item[1] + ) + else: # it gives it the mat ID + for i in range(len(mcdc["materials"])): + legend_hands[i] = mpatches.Patch(color=colors[i], label="{}".format(i)) + # Set reference axis for axis in ["x", "y", "z"]: if axis not in vis_type: @@ -2498,6 +2735,8 @@ def visualize( particle["g"] = 0 particle["E"] = 1e6 + frames = [] + for t in time: # Set time particle["t"] = t @@ -2516,6 +2755,7 @@ def visualize( for j in range(pixels[1]): particle[second_key] = second_midpoint[j] + # print( particle['x'], particle['y'], particle['z']) # Get material particle["cell_ID"] = -1 particle["material_ID"] = -1 @@ -2529,8 +2769,16 @@ def visualize( plt.xlabel(first_key + " [cm]") plt.ylabel(second_key + " [cm]") plt.title(reference_key + " = %.2f cm" % reference + ", time = %.2f s" % t) + if legend: + plt.legend(handles=legend_hands, bbox_to_anchor=(1.04, 1), loc="upper left") if save_as is not None: plt.savefig(save_as + "_%.2f.png" % t) plt.clf() + if movie: + image = imageio.v2.imread(save_as + "_%.2f.png" % t) + frames.append(image) else: plt.show() + + if movie: + imageio.mimsave(save_as + ".gif", frames, fps=5) From 5ef559b986e2c389b441bd153a91e0055dd51277 Mon Sep 17 00:00:00 2001 From: "Joanna Piper Morgan (jonsey)" Date: Thu, 24 Apr 2025 16:11:08 -0700 Subject: [PATCH 10/17] adding notes for delta tracking --- mcdc/print_.py | 21 +++++++++++++++++++-- 1 file changed, 19 insertions(+), 2 deletions(-) diff --git a/mcdc/print_.py b/mcdc/print_.py index 9ea000cb..b94b3844 100644 --- a/mcdc/print_.py +++ b/mcdc/print_.py @@ -26,7 +26,7 @@ def print_warning(msg): sys.stdout.flush() -def print_banner(mcdc): +def print_banner(mcdc, args): size = MPI.COMM_WORLD.Get_size() if master: banner = ( @@ -47,6 +47,11 @@ def print_banner(mcdc): banner += " Mode | Python\n" else: banner += " Mode | Numba\n" + if args.target == "cpu": + banner += " Target | CPU\n" + elif args.target == "gpu": + banner += " Target | GPU\n" + if mcdc["technique"]["iQMC"]: banner += " Algorithm | iQMC\n" if mcdc["setting"]["mode_eigenvalue"]: @@ -54,10 +59,22 @@ def print_banner(mcdc): else: solver = mcdc["technique"]["iqmc"]["fixed_source_solver"] banner += " Solver | " + solver + "\n" + elif ( + mcdc["technique"]["delta_tracking"] + and not mcdc["technique"]["collision_estimator"] + ): + banner += " Algorithm | History-based\n" + banner += " | hybrid delta tracking\n" + elif ( + mcdc["technique"]["delta_tracking"] + and mcdc["technique"]["collision_estimator"] + ): + banner += " Algorithm | History-based\n" + banner += " | delta tracking (collision_estimator)\n" else: banner += " Algorithm | History-based\n" banner += " MPI Processes | %i\n" % size - banner += " OpenMP Threads | 1" + # banner += " OpenMP Threads | 1" print(banner) sys.stdout.flush() From 088800ef086b464a1572153454c2abf79c34ae04 Mon Sep 17 00:00:00 2001 From: "Joanna Piper Morgan (jonsey)" Date: Thu, 24 Apr 2025 16:11:49 -0700 Subject: [PATCH 11/17] adding a lookup to boundary operaiton and removing lost particle reporting for delta tracking --- mcdc/src/geometry.py | 39 +++++++++++++++++++++++++++++++++++++-- 1 file changed, 37 insertions(+), 2 deletions(-) diff --git a/mcdc/src/geometry.py b/mcdc/src/geometry.py index b8b20899..cea3ab1e 100644 --- a/mcdc/src/geometry.py +++ b/mcdc/src/geometry.py @@ -254,9 +254,14 @@ def locate_particle(particle_container, mcdc): particle["uy"] = uy_global particle["uz"] = uz_global - # Report lost particle + # Report lost particle if not delta tracking if particle_is_lost: - report_lost(particle_container) + # Cheaper then keeping track of vac bounds + # reflecting boundaries are still tracked + if particle["delta_tracking"]: + particle["alive"] = False + else: + report_lost(particle_container) return not particle_is_lost @@ -451,6 +456,36 @@ def distance_to_nearest_surface(particle_container, cell, mcdc): return distance, surface_ID +@njit +def distance_to_boundary(particle_container, mcdc): + """ + Determine the nearest cell surface and the distance to it + """ + + particle = particle_container[0] + distance = INF + surface_ID = -1 + + # Particle parameters + speed = physics.get_speed(particle_container, mcdc) + N_surface = len(mcdc["boundary_surface_IDs"]) # array of ints + + for surf in range(N_surface): + candidate_surface_ID = mcdc["boundary_surface_IDs"][surf] + surface = mcdc["surfaces"][candidate_surface_ID] + d = surface_.get_distance(particle_container, speed, surface) + if d < distance: + distance = d + surface_ID = surface["ID"] + + # if it's an int it never hits a surface! + if distance is not INF: + particle["surface_ID"] = surface_ID + particle["event"] = EVENT_SURFACE_CROSSING + + return distance + + # ====================================================================================== # Miscellanies # ====================================================================================== From 83f196d1f69f522fc03b59adbc09103011dc991e Mon Sep 17 00:00:00 2001 From: "Joanna Piper Morgan (jonsey)" Date: Thu, 24 Apr 2025 16:12:26 -0700 Subject: [PATCH 12/17] adding imageio for movie making in the visulizer --- pyproject.toml | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/pyproject.toml b/pyproject.toml index f3881efd..bb6fa0f8 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -59,7 +59,8 @@ dependencies = [ "black", "sympy", "pre-commit", - "cvxpy" + "cvxpy", + "imageio", ] [project.optional-dependencies] docs = ["sphinx==7.2.6", "furo", "sphinx_toolbox"] From 127d80f8f671145f35bf4ab785fc991dc82b1568 Mon Sep 17 00:00:00 2001 From: "Joanna Piper Morgan (jonsey)" Date: Thu, 24 Apr 2025 16:16:39 -0700 Subject: [PATCH 13/17] adding delta tracking to input def --- docs/source/pythonapi/index.rst | 1 + 1 file changed, 1 insertion(+) diff --git a/docs/source/pythonapi/index.rst b/docs/source/pythonapi/index.rst index 8cb0b30b..7403cd61 100644 --- a/docs/source/pythonapi/index.rst +++ b/docs/source/pythonapi/index.rst @@ -63,6 +63,7 @@ Defining techniques mcdc.weighted_emission mcdc.weight_roulette mcdc.weight_window + mcdc.delta_tracking From 61768ca59b0bd7196b0429c4cf2a707990031261 Mon Sep 17 00:00:00 2001 From: "Joanna Piper Morgan (jonsey)" Date: Sat, 26 Apr 2025 23:05:04 -0700 Subject: [PATCH 14/17] removing false interpretation of delta tracking surface tracking switch condition --- mcdc/kernel.py | 7 +------ 1 file changed, 1 insertion(+), 6 deletions(-) diff --git a/mcdc/kernel.py b/mcdc/kernel.py index 06e645ce..49fa0947 100644 --- a/mcdc/kernel.py +++ b/mcdc/kernel.py @@ -2856,10 +2856,6 @@ def move_to_event(P_arr, data, mcdc): # don't delta track if in excluded material if P["material_ID"] == mcdc["technique"]["dt_material_exclusion"]: P["delta_tracking"] = False - - # but complete rejection sampleing process before switching - if P["event"] & EVENT_PHANTOM_COLLISION: - P["delta_tracking"] = True else: P["delta_tracking"] = False @@ -2998,8 +2994,7 @@ def move_to_event(P_arr, data, mcdc): # Mesh tallies # use a collision estimator for currently undefined tallies for tally in mcdc["mesh_tallies"]: - if tally is not SCORE_FLUX: - score_mesh_collision_tally(P_arr, tally, data, mcdc) + score_mesh_collision_tally(P_arr, tally, data, mcdc) @njit From faf5029700a64fbdb873b587fb27fb5667815d0b Mon Sep 17 00:00:00 2001 From: Joanna Piper Morgan Date: Sun, 4 May 2025 13:58:09 -0700 Subject: [PATCH 15/17] delta tracking never reports lost particles --- mcdc/src/geometry.py | 8 +++++++- 1 file changed, 7 insertions(+), 1 deletion(-) diff --git a/mcdc/src/geometry.py b/mcdc/src/geometry.py index cea3ab1e..1985e1e4 100644 --- a/mcdc/src/geometry.py +++ b/mcdc/src/geometry.py @@ -152,7 +152,13 @@ def inspect_geometry(particle_container, mcdc): # Report lost particle if event == EVENT_LOST: - report_lost(particle_container) + # Cheaper then keeping track of vac bounds + # reflecting boundaries are still tracked + # needed in this function for hybrid methods + if particle["delta_tracking"]: + particle["alive"] = False + else: + report_lost(particle_container) # Assign particle event particle["event"] = event From b30e25a4bc18e94eba755bf9b96250fd721db3e4 Mon Sep 17 00:00:00 2001 From: Joanna Piper Morgan Date: Mon, 5 May 2025 14:57:12 -0700 Subject: [PATCH 16/17] switching particle killings back to technique based --- mcdc/src/geometry.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/mcdc/src/geometry.py b/mcdc/src/geometry.py index 1985e1e4..3cc4ee1d 100644 --- a/mcdc/src/geometry.py +++ b/mcdc/src/geometry.py @@ -155,7 +155,7 @@ def inspect_geometry(particle_container, mcdc): # Cheaper then keeping track of vac bounds # reflecting boundaries are still tracked # needed in this function for hybrid methods - if particle["delta_tracking"]: + if mcdc["technique"]["delta_tracking"]: particle["alive"] = False else: report_lost(particle_container) @@ -264,7 +264,7 @@ def locate_particle(particle_container, mcdc): if particle_is_lost: # Cheaper then keeping track of vac bounds # reflecting boundaries are still tracked - if particle["delta_tracking"]: + if mcdc["technique"]["delta_tracking"]: particle["alive"] = False else: report_lost(particle_container) From ce74a458d5f0e96cd652f9e2cef1583085952dc7 Mon Sep 17 00:00:00 2001 From: "Joanna Piper Morgan (jonsey)" Date: Wed, 7 May 2025 13:57:43 -0700 Subject: [PATCH 17/17] fixing collision estimator for other reaction desnities --- mcdc/kernel.py | 8 ++------ 1 file changed, 2 insertions(+), 6 deletions(-) diff --git a/mcdc/kernel.py b/mcdc/kernel.py index 49fa0947..bbf6f4c6 100644 --- a/mcdc/kernel.py +++ b/mcdc/kernel.py @@ -1851,7 +1851,7 @@ def mesh_get_energy_index(P_arr, mesh, mode_MG): @njit def score_mesh_collision_tally(P_arr, tally, data, mcdc): - """Scoring tallies with a collision estimator,""" + """Scoring tallies with a collision estimator""" P = P_arr[0] tally_bin = data[TALLY] material = mcdc["materials"][P["material_ID"]] @@ -1915,10 +1915,6 @@ def score_mesh_collision_tally(P_arr, tally, data, mcdc): # score = flux * mu * mu * (x - (mesh["x"][ix + 1] + mesh["x"][ix]) / 2) adapt.global_add(tally_bin, (TALLY_SCORE, idx + i), round(score)) - SigmaT = get_MacroXS(XS_TOTAL, material, P_arr, mcdc) - score = 1.0 / SigmaT - adapt.global_add(tally_bin, (TALLY_SCORE, idx + i), round(score)) - return @@ -2990,7 +2986,7 @@ def move_to_event(P_arr, data, mcdc): rejection_sample(P_arr, mcdc) if mcdc["technique"]["collision_estimator"]: - if P["event"] == EVENT_COLLISION and mcdc["cycle_active"]: + if (P["event"] == EVENT_COLLISION or P["event"]== EVENT_PHANTOM_COLLISION) and mcdc["cycle_active"]: # Mesh tallies # use a collision estimator for currently undefined tallies for tally in mcdc["mesh_tallies"]: