diff --git a/docs/source/pythonapi/index.rst b/docs/source/pythonapi/index.rst index 8cb0b30b1..7403cd612 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 diff --git a/mcdc/__init__.py b/mcdc/__init__.py index 6851a6f44..ec7a9eaf1 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, diff --git a/mcdc/config.py b/mcdc/config.py index 3b80f3d04..4877efad6 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__)) diff --git a/mcdc/constant.py b/mcdc/constant.py index 1667f65cd..ec96afe45 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 diff --git a/mcdc/global_.py b/mcdc/global_.py index 64895cca7..f062df314 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 = { diff --git a/mcdc/input_.py b/mcdc/input_.py index fa29af356..080ef9778 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. diff --git a/mcdc/kernel.py b/mcdc/kernel.py index 436a7d9fa..fbff82bea 100644 --- a/mcdc/kernel.py +++ b/mcdc/kernel.py @@ -1849,6 +1849,75 @@ 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)) + + return + + @njit def score_mesh_tally(P_arr, distance, tally, data_tally, mcdc): P = P_arr[0] @@ -2757,11 +2826,29 @@ def move_to_event(P_arr, data_tally, 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 + 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 @@ -2775,12 +2862,19 @@ def move_to_event(P_arr, data_tally, 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 # ================================================================================== @@ -2840,40 +2934,64 @@ def move_to_event(P_arr, data_tally, 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_tally, 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_tally, mcdc) + # ========================================================================= + # Delta Tracking + # ========================================================================= - # CS tallies - for tally in mcdc["cs_tallies"]: - score_cs_tally(P_arr, distance, tally, data_tally, 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 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"]: + 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: @@ -2882,9 +3000,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 # ============================================================================= @@ -3699,6 +3837,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 # ============================================================================= @@ -3746,6 +3924,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] diff --git a/mcdc/loop.py b/mcdc/loop.py index d2a8105f2..973d12e4b 100644 --- a/mcdc/loop.py +++ b/mcdc/loop.py @@ -562,7 +562,6 @@ def gpu_loop_source(seed, data, mcdc): def loop_particle(P_arr, data_tally, prog): P = P_arr[0] mcdc = adapt.mcdc_global(prog) - while P["alive"]: step_particle(P_arr, data_tally, prog) @@ -577,6 +576,7 @@ def step_particle(P_arr, data_tally, prog): # Execute events if P["event"] == EVENT_LOST: + P["alive"] = False return # Collision @@ -604,6 +604,10 @@ def step_particle(P_arr, data_tally, 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_tally, prog) diff --git a/mcdc/main.py b/mcdc/main.py index 63558206e..5f56e517d 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() @@ -69,7 +70,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"]: @@ -643,7 +644,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] @@ -693,8 +696,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"]: @@ -1220,6 +1227,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) @@ -1442,6 +1453,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 # ========================================================================= @@ -2412,7 +2611,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, ): """ @@ -2434,9 +2636,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] @@ -2452,6 +2675,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: @@ -2500,6 +2737,8 @@ def visualize( particle["g"] = 0 particle["E"] = 1e6 + frames = [] + for t in time: # Set time particle["t"] = t @@ -2518,6 +2757,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 @@ -2531,8 +2771,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) diff --git a/mcdc/print_.py b/mcdc/print_.py index 9ea000cb4..b94b38442 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() diff --git a/mcdc/src/geometry.py b/mcdc/src/geometry.py index b8b208996..3cc4ee1d1 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 mcdc["technique"]["delta_tracking"]: + particle["alive"] = False + else: + report_lost(particle_container) # Assign particle event particle["event"] = event @@ -254,9 +260,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 mcdc["technique"]["delta_tracking"]: + particle["alive"] = False + else: + report_lost(particle_container) return not particle_is_lost @@ -451,6 +462,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 # ====================================================================================== diff --git a/mcdc/type_.py b/mcdc/type_.py index b86bc1696..240f6dbbe 100644 --- a/mcdc/type_.py +++ b/mcdc/type_.py @@ -220,6 +220,7 @@ def make_type_particle(input_deck): ("fresh", bool_), ("event", int64), ("rng_seed", uint64), + ("delta_tracking", bool_), ] # Get modes @@ -1099,6 +1100,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 ] # ========================================================================= @@ -1286,10 +1291,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 @@ -1448,6 +1481,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]) @@ -1514,6 +1555,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,)), diff --git a/pyproject.toml b/pyproject.toml index 335257fda..61ab26a6b 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"]