From 86560e70c79c0be7d9a547c684db5be134b30190 Mon Sep 17 00:00:00 2001 From: Jackson Morgan <67020956+jpmorgan98@users.noreply.github.com> Date: Tue, 7 Feb 2023 12:08:09 -0800 Subject: [PATCH 1/9] Update README.md --- README.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/README.md b/README.md index d945675e..47fdbc0d 100644 --- a/README.md +++ b/README.md @@ -27,7 +27,7 @@ on the [AZURV1 benchmark](https://inis.iaea.org/search/search.aspx?orig_q=RN:410 ```python import numpy as np -import mcdc +import mcdc test # ============================================================================= # Set model From 235fd9983eef81e685efa08932ce49a4991b5f24 Mon Sep 17 00:00:00 2001 From: Jackson Morgan <67020956+jpmorgan98@users.noreply.github.com> Date: Wed, 1 Mar 2023 17:11:02 -0800 Subject: [PATCH 2/9] Update README.md --- README.md | 2 ++ 1 file changed, 2 insertions(+) diff --git a/README.md b/README.md index 98e511f5..60d0aea5 100644 --- a/README.md +++ b/README.md @@ -1,5 +1,7 @@ # MC/DC: Monte Carlo Dynamic Code +## test alteration do not approve the PR! + ![mcdc_logo v1](https://user-images.githubusercontent.com/26186244/173467190-74d9b09a-ef7d-4f0e-8bdf-4a076de7c43c.svg) [![License](https://img.shields.io/badge/License-BSD_3--Clause-blue.svg)](https://opensource.org/licenses/BSD-3-Clause) From 0db880b8201f6f502cdafe5776c268c64b73bb96 Mon Sep 17 00:00:00 2001 From: Joanna Piper Morgan Date: Thu, 9 Mar 2023 12:20:10 -0800 Subject: [PATCH 3/9] added stuff to readme --- README.md | 3 +++ 1 file changed, 3 insertions(+) diff --git a/README.md b/README.md index 60d0aea5..f26b83f2 100644 --- a/README.md +++ b/README.md @@ -12,6 +12,9 @@ MC/DC is a performant, scalable, and machine-portable Python-based Monte Carlo neutron transport software currently developed in the Center for Exascale Monte Carlo Neutron Transport ([CEMeNT](https://cement-psaap.github.io/)). + +This is a test blah blah blah + *Please Note that this project is in the early stages of devlopment (not even an alpha).* That being said enjoy! ## Installation From 45f8d0044c6de5b8177c22612af31d72663bcd89 Mon Sep 17 00:00:00 2001 From: Jackson Morgan <67020956+jpmorgan98@users.noreply.github.com> Date: Tue, 23 May 2023 11:48:40 -0700 Subject: [PATCH 4/9] Update README.md --- README.md | 2 -- 1 file changed, 2 deletions(-) diff --git a/README.md b/README.md index 2176c954..a814fcde 100644 --- a/README.md +++ b/README.md @@ -1,7 +1,5 @@ # MC/DC: Monte Carlo Dynamic Code -## test alteration do not approve the PR! - ![mcdc_logo v1](https://user-images.githubusercontent.com/26186244/173467190-74d9b09a-ef7d-4f0e-8bdf-4a076de7c43c.svg) [![License](https://img.shields.io/badge/License-BSD_3--Clause-blue.svg)](https://opensource.org/licenses/BSD-3-Clause) From 5728ca1ccf4c9b6524f417976682e44c68ba6935 Mon Sep 17 00:00:00 2001 From: "Joanna Piper Morgan (jonsey)" Date: Thu, 21 Dec 2023 17:25:54 -0800 Subject: [PATCH 5/9] Mergin visulization functionality previously implemented by Rhoan. Passing tests. Adding new file to /mcdc visulizer.py --- .github/workflows/tests.yml | 2 +- examples/godiva-mockup-visulization/input.py | 96 +++++ .../godiva-mockup-visulization/process.py | 41 ++ mcdc/__init__.py | 1 + mcdc/card.py | 3 + mcdc/input_.py | 9 + mcdc/visualizer.py | 403 ++++++++++++++++++ 7 files changed, 554 insertions(+), 1 deletion(-) create mode 100644 examples/godiva-mockup-visulization/input.py create mode 100644 examples/godiva-mockup-visulization/process.py create mode 100644 mcdc/visualizer.py diff --git a/.github/workflows/tests.yml b/.github/workflows/tests.yml index e224c51f..a1447582 100644 --- a/.github/workflows/tests.yml +++ b/.github/workflows/tests.yml @@ -19,7 +19,7 @@ jobs: - name: Install dependencies run: | #sudo apt-get install libopenmpi-dev - pip install numpy numba h5py matplotlib scipy pytest colorama mpi4py + pip install numpy numba h5py matplotlib scipy pytest colorama mpi4py ngsolve distinctipy pip list pip install -e . - name: Patch Numba diff --git a/examples/godiva-mockup-visulization/input.py b/examples/godiva-mockup-visulization/input.py new file mode 100644 index 00000000..45148935 --- /dev/null +++ b/examples/godiva-mockup-visulization/input.py @@ -0,0 +1,96 @@ +import numpy as np +import mcdc, h5py + +# ============================================================================= +# Materials +# ============================================================================= + +m_abs = mcdc.material(capture=np.array([1e5]), speed=np.array([1e3]), name="water") +m_void = mcdc.material( + capture=np.array([5e-5]), + scatter=np.array([[5e-5]]), + speed=np.array([1e3]), + name="source", +) + +# ============================================================================= +# Set surfaces +# ============================================================================= + +# For cube boundaries +cube_x0 = mcdc.surface("plane-x", x=-22.0, bc="vacuum") +cube_x1 = mcdc.surface("plane-x", x=22.0, bc="vacuum") +cube_y0 = mcdc.surface("plane-y", y=-12.0, bc="vacuum") +cube_y1 = mcdc.surface("plane-y", y=12.0, bc="vacuum") +cube_z0 = mcdc.surface("plane-z", z=-12.0, bc="vacuum") +cube_z1 = mcdc.surface("plane-z", z=12.0, bc="vacuum") + +# For the 3-part hollow sphere +sp_left = mcdc.surface("sphere", center=[-2.0, 0.0, 0.0], radius=6.0) +sp_center = mcdc.surface("sphere", center=[0.0, 0.0, 0.0], radius=6.0) +sp_right = mcdc.surface("sphere", center=[2.0, 0.0, 0.0], radius=6.0) +pl_x0 = mcdc.surface("plane-x", x=-3.5) +pl_x1 = mcdc.surface("plane-x", x=-1.5) +pl_x2 = mcdc.surface("plane-x", x=1.5) +pl_x3 = mcdc.surface("plane-x", x=3.5) + +# For the moving rod +cy = mcdc.surface("cylinder-x", center=[0.0, 0.0], radius=0.5) +pl_rod0 = mcdc.surface("plane-x", x=[-22.0, 22.0 - 12.0], t=[0.0, 5.0]) +pl_rod1 = mcdc.surface("plane-x", x=[-22.0 + 12.0, 22.0], t=[0.0, 5.0]) + +# ============================================================================= +# Set cells +# ============================================================================= + +# Moving rod +mcdc.cell([-cy, +pl_rod0, -pl_rod1], m_void) + +# 3-part hollow shpere +mcdc.cell([-sp_left, -pl_x0, +cy], m_void) +mcdc.cell([-sp_center, +pl_x1, -pl_x2, +cy], m_void) +mcdc.cell([-sp_right, +pl_x3, +cy], m_void) + +# Surrounding water +# Left of rod +mcdc.cell([-cy, +cube_x0, -pl_rod0], m_abs) +# Right of rod +mcdc.cell([-cy, +pl_rod1, -cube_x1], m_abs) +# The rest +mcdc.cell( + [+cy, +sp_left, +cube_x0, -pl_x0, +cube_y0, -cube_y1, +cube_z0, -cube_z1], m_abs +) +mcdc.cell([+cy, +pl_x0, -pl_x1, +cube_y0, -cube_y1, +cube_z0, -cube_z1], m_abs) +mcdc.cell([+sp_center, +pl_x1, -pl_x2, +cube_y0, -cube_y1, +cube_z0, -cube_z1], m_abs) +mcdc.cell([+cy, +pl_x2, -pl_x3, +cube_y0, -cube_y1, +cube_z0, -cube_z1], m_abs) +mcdc.cell( + [+cy, +sp_right, +pl_x3, -cube_x1, +cube_y0, -cube_y1, +cube_z0, -cube_z1], m_abs +) + +# ============================================================================= +# Set source +# ============================================================================= + +mcdc.source(x=[-22.0, 22.0], time=[0.0, 5.0], isotropic=True) + +mcdc.visualize(start_time=0, end_time=5) +# ============================================================================= +# Set tally, setting, and run mcdc +# ============================================================================= + +# Tally: z-integrated flux (X-Y section view) + +""" +mcdc.tally( + scores=["flux"], + x=np.linspace(-22.0, 22.0, 84+1), + y=np.linspace(-12.0, 12.0, 24+1), + t=np.linspace(0.0, 5.0, 50+1), +) + +# Setting +mcdc.setting(N_particle=1e6) + +# Run +mcdc.run() +""" diff --git a/examples/godiva-mockup-visulization/process.py b/examples/godiva-mockup-visulization/process.py new file mode 100644 index 00000000..78280fe9 --- /dev/null +++ b/examples/godiva-mockup-visulization/process.py @@ -0,0 +1,41 @@ +import numpy as np +import matplotlib.pyplot as plt +from matplotlib.colors import LogNorm +import h5py +import matplotlib.animation as animation + + +# ============================================================================= +# Plot results +# ============================================================================= + +# Results +with h5py.File("output.h5", "r") as f: + x = f["tally/grid/x"][:] + x_mid = 0.5 * (x[:-1] + x[1:]) + y = f["tally/grid/y"][:] + y_mid = 0.5 * (y[:-1] + y[1:]) + t = f["tally/grid/t"][:] + t_mid = 0.5 * (t[:-1] + t[1:]) + X, Y = np.meshgrid(y, x) + + phi = f["tally/flux/mean"][:] + phi_sd = f["tally/flux/sdev"][:] + +fig, ax = plt.subplots() +cax = ax.pcolormesh(X, Y, phi[0], vmin=phi[0].min(), vmax=phi[0].max()) +text = ax.text(0.02, 1.02, "", transform=ax.transAxes) +ax.set_aspect("equal", "box") +ax.set_xlabel("$y$ [cm]") +ax.set_ylabel("$x$ [cm]") +ax.set_aspect("equal") + + +def animate(i): + cax.set_array(phi[i]) + cax.set_clim(phi[i].min(), phi[i].max()) + text.set_text(r"$t \in [%.1f,%.1f]$ s" % (t[i], t[i + 1])) + + +anim = animation.FuncAnimation(fig, animate, interval=10, frames=len(t) - 1) +plt.show() diff --git a/mcdc/__init__.py b/mcdc/__init__.py index bc43b278..09595dfb 100644 --- a/mcdc/__init__.py +++ b/mcdc/__init__.py @@ -25,3 +25,4 @@ reset_cards, ) from mcdc.main import run, prepare +from mcdc.visualizer import visualize diff --git a/mcdc/card.py b/mcdc/card.py index ca69c398..05b26a08 100644 --- a/mcdc/card.py +++ b/mcdc/card.py @@ -212,6 +212,7 @@ def make_card_material(N_nuclide, G, J): card["nu_f"] = np.zeros(G) card["chi_s"] = np.zeros([G, G]) card["chi_p"] = np.zeros([G, G]) + card["name"] = None card["sensitivity"] = False card["uq"] = False return card @@ -241,6 +242,7 @@ def make_card_surface(): card["nz"] = 0.0 card["sensitivity"] = False card["sensitivity_ID"] = 0 + card["type"] = " " card["dsm_Np"] = 1.0 return card @@ -253,6 +255,7 @@ def make_card_cell(N_surface): card["surface_IDs"] = np.zeros(N_surface, dtype=int) card["positive_flags"] = np.zeros(N_surface, dtype=bool) card["material_ID"] = 0 + card["material_name"] = None card["lattice"] = False card["lattice_ID"] = 0 card["lattice_center"] = np.array([0.0, 0.0, 0.0]) diff --git a/mcdc/input_.py b/mcdc/input_.py index 8ab1b350..150602e9 100644 --- a/mcdc/input_.py +++ b/mcdc/input_.py @@ -229,6 +229,7 @@ def material( chi_d=None, speed=None, decay=None, + name="P", sensitivity=False, dsm_Np=1.0, ): @@ -306,6 +307,11 @@ def material( card = make_card_material(N_nuclide, G, J) card["ID"] = len(mcdc.input_deck.materials) + if name is not None: + card["name"] = name + else: + card["name"] = card["ID"] + # Calculate basic XS and determine sensitivity flag for i in range(N_nuclide): nuc = nuclides[i][0] @@ -476,6 +482,8 @@ def surface(type_, bc="interface", sensitivity=False, dsm_Np=1.0, **kw): # Axx + Byy + Czz + Dxy + Exz + Fyz + Gx + Hy + Iz + J(t) = 0 # J(t) = J0_i + J1_i*t for t in [t_{i-1}, t_i), t_0 = 0 + card["type"] = type_ + # Set up surface attributes if type_ == "plane-x": check_requirement("surface plane-x", kw, ["x"]) @@ -628,6 +636,7 @@ def cell(surfaces_flags, fill, lattice_center=None): # Material cell else: card["material_ID"] = fill["ID"] + card["material_name"] = fill["name"] # Push card mcdc.input_deck.cells.append(card) diff --git a/mcdc/visualizer.py b/mcdc/visualizer.py new file mode 100644 index 00000000..656c8bb2 --- /dev/null +++ b/mcdc/visualizer.py @@ -0,0 +1,403 @@ +from netgen.meshing import * +from netgen.csg import * +from ngsolve import Draw, Redraw # just for visualization +from tkinter import * # Tkinter is used to create the window for the time slider and color key +import distinctipy # creates unlimited visually distinct colors for visualization +import math + +# Get input_card and set global variables as "mcdc_" + +import mcdc.global_ as mcdc_ + +input_card = mcdc_.input_deck + + +# get a point on the plane based on the current time in the system +# start and end times are zero by default +# called with a plane is visualized in create_cell_geometry() +def get_plane_current_position(surface, current_time, start_time, end_time): + if len(surface["t"]) > 2: # check if shape moves + # establish reference points + start_time = start_time # default start time is zero + end_time = end_time + start_position = -surface["J"][0][0] + end_position = -surface["J"][1][0] + + current_position = start_position + ( + ((end_position - start_position) / (end_time - start_time)) * current_time + ) + + else: + current_position = -surface["J"][0][0] + + return current_position + + +# create the CSG geometry for a cell +# called by draw_geometry() +def create_cell_geometry(cell, current_time, surface_list, start_time, end_time): + cell_shape_list = [] # list of shapes that make up the cell + for i in range(0, len(cell["surface_IDs"])): + surface_ID = cell["surface_IDs"][i] + + ## Check type of shape + + if surface_list[surface_ID]["type"] == "plane-x": + # get reference point from the surface card + point = ( + ( + get_plane_current_position( + surface_list[surface_ID], + current_time, + start_time=start_time, + end_time=end_time, + ) + ), + 0, + 0, + ) + + # get normal vector from the surface card + if cell["positive_flags"][i]: + vector = ( + -surface_list[surface_ID]["G"], + surface_list[surface_ID]["H"], + surface_list[surface_ID]["I"], + ) + else: + vector = ( + surface_list[surface_ID]["G"], + surface_list[surface_ID]["H"], + surface_list[surface_ID]["I"], + ) + + # planes have to be intersected to achive the wanted visualization in ngsolve + cell_shape_list.append([Plane(Pnt(point), Vec(vector)), "intersect"]) + + elif surface_list[surface_ID]["type"] == "plane-y": + # get reference point from the surface card + point = ( + 0, + ( + get_plane_current_position( + surface_list[surface_ID], + current_time, + start_time=start_time, + end_time=end_time, + ) + ), + 0, + ) + + # get normal vector from the surface card + if cell["positive_flags"][i]: + vector = ( + surface_list[surface_ID]["G"], + -surface_list[surface_ID]["H"], + surface_list[surface_ID]["I"], + ) + else: + vector = ( + surface_list[surface_ID]["G"], + surface_list[surface_ID]["H"], + surface_list[surface_ID]["I"], + ) + + # planes have to be intersected to achive the wanted visualization in ngsolve + cell_shape_list.append( + [Plane(Pnt(point), Vec(vector)).col([1, 0, 0]), "intersect"] + ) + + elif surface_list[surface_ID]["type"] == "plane-z": + # get reference point from the surface card + point = ( + 0, + 0, + ( + get_plane_current_position( + surface_list[surface_ID], + current_time, + start_time=start_time, + end_time=end_time, + ) + ), + ) + + # get normal vector from the surface card + if cell["positive_flags"][i]: + vector = ( + surface_list[surface_ID]["G"], + surface_list[surface_ID]["H"], + -surface_list[surface_ID]["I"], + ) + else: + vector = ( + surface_list[surface_ID]["G"], + surface_list[surface_ID]["H"], + surface_list[surface_ID]["I"], + ) + + # planes have to be intersected to achive the wanted visualization in ngsolve + cell_shape_list.append( + [Plane(Pnt(point), Vec(vector)).col([1, 0, 0]), "intersect"] + ) + + elif surface_list[surface_ID]["type"] == "sphere": + # Get the center point from the surface card + x = surface_list[surface_ID]["G"] / -2 + y = surface_list[surface_ID]["H"] / -2 + z = surface_list[surface_ID]["I"] / -2 + point = (x, y, z) + + # get radius from the surface card + radius = float( + math.sqrt( + abs(surface_list[surface_ID]["J"][0, 0] - x**2 - y**2 - z**2) + ) + ) + + # Add or subtract the sphere based on the CSG input + if cell["positive_flags"][i]: + cell_shape_list.append([Sphere(Pnt(point), radius), "subtract"]) + + else: + cell_shape_list.append([Sphere(Pnt(point), radius), "add"]) + + elif surface_list[surface_ID]["type"] == "cylinder-x": + # the y and z points are used to define the central axis of the cylinder + y = surface_list[surface_ID]["H"] / -2 + z = surface_list[surface_ID]["I"] / -2 + + # radius of the circle formed by a cross section parallel to the yz-plane + radius = float( + math.sqrt(abs(surface_list[surface_ID]["J"][0, 0] - z**2 - y**2)) + ) + + # Add or subtract the sphere based on the CSG input + # in the NGSolve CSG cylinders are infinite along their central axis, + # by default we make the central axis -100<=x<=100 + # if a larger central axis is needed change below + if cell["positive_flags"][i]: + cell_shape_list.append( + [Cylinder(Pnt(-100, y, z), Pnt(100, y, z), radius), "subtract"] + ) + else: + cell_shape_list.append( + [Cylinder(Pnt(-100, y, z), Pnt(100, y, z), radius), "add"] + ) + elif surface_list[surface_ID]["type"] == "cylinder-y": + # the x and z points are used to define the central axis of the cylinder + x = surface_list[surface_ID]["G"] / -2 + z = surface_list[surface_ID]["I"] / -2 + + # radius of the circle formed by a cross section parallel to the xz-plane + radius = float( + math.sqrt(abs(surface_list[surface_ID]["J"][0, 0] - z**2 - y**2)) + ) + + # Add or subtract the sphere based on the CSG input + # in the NGSolve CSG cylinders are infinite along their central axis, + # by default we make the central axis -100<=y<=100 + # if a larger central axis is needed change below + if cell["positive_flags"][i]: + cell_shape_list.append( + [Cylinder(Pnt(x, -100, z), Pnt(x, 100, z), radius), "subtract"] + ) + else: + cell_shape_list.append( + [Cylinder(Pnt(x, -100, z), Pnt(x, 100, z), radius), "add"] + ) + elif surface_list[surface_ID]["type"] == "cylinder-z": + # the x and y points are used to define the central axis of the cylinder + x = surface_list[surface_ID]["G"] / -2 + y = surface_list[surface_ID]["H"] / -2 + + # radius of the circle formed by a cross section parallel to the xy-plane + radius = float( + math.sqrt(abs(surface_list[surface_ID]["J"][0, 0] - z**2 - y**2)) + ) + + # Add or subtract the sphere based on the CSG input + # in the NGSolve CSG cylinders are infinite along their central axis, + # by default we make the central axis -100<=z<=100 + # if a larger central axis is needed change below + if cell["positive_flags"][i]: + cell_shape_list.append( + [Cylinder(Pnt(x, y, -100), Pnt(x, y, 100), radius), "subtract"] + ) + else: + cell_shape_list.append( + [Cylinder(Pnt(x, y, -100), Pnt(x, y, 100), radius), "add"] + ) + + # combine the shapes for each cell + # for subtraction to work properly, it must be done in the end + for i in range(0, len(cell_shape_list)): + if not (cell_shape_list[i][1] == "subtract"): + if "cell_geometry" not in locals(): + cell_geometry = cell_shape_list[i][0] + elif cell_shape_list[i][1] == "intersect": + cell_geometry = cell_geometry * cell_shape_list[i][0] + + elif cell_shape_list[i][1] == "add": + cell_geometry = cell_geometry + cell_shape_list[i][0] + + for i in range(0, len(cell_shape_list)): + if cell_shape_list[i][1] == "subtract": + cell_geometry = cell_geometry - cell_shape_list[i][0] + + return cell_geometry # return the finished CSG geometry + + +# visualizes the model at a specified time (current_time, type float) +# called by visualize() +def draw_Geometry(current_time, start_time, end_time): + # create lists that contain all cells and surfaces + surface_list = input_card.surfaces + cell_list = input_card.cells + # surface_list = input_card.surfaces + + # make water blue and the source green + water_rgb = [0, 0, 1] + source_rgb = [0, 1, 0] + + geo = CSGeometry() # create the ngsolve geometry object + + # list of materials that need colors to be generated (ie not water or the source) + material_colors_to_generate = [] + + # find the materials that need colors generated and add them to material_colors_to_generate + for cell in cell_list: + cell_material_name = cell["material_name"] + if ( + (cell_material_name not in material_colors_to_generate) + and (cell_material_name != "water") + and (cell_material_name != "source") + ): + material_colors_to_generate.append(cell_material_name) + + # colors that should not be generated (ie, taken by preset materials or which are visually unappealing) + # These colors are rgb values, more can be added by extending the list + input_colors = [ + (water_rgb[0], water_rgb[1], water_rgb[2]), # water - blue + (source_rgb[0], source_rgb[1], source_rgb[2]), # source - green + (1, 1, 1), # white + (0, 0, 0), # black + ] + # create n number of distinct colors where n is the number of materials in material_colors_to_generate + distinct_colors = distinctipy.get_colors( + len(material_colors_to_generate), input_colors + ) + + # This list will later be passed to create_color_key + # contains lists of format [rgb value, material name] + color_key_list = [] + materials_added_to_color_key = [] + + # cycle through the cells in the model + for cell_index in range(0, len(cell_list)): + cell = cell_list[cell_index] + + # create the geometry for the cell + cell_geometry = create_cell_geometry( + cell, + current_time=current_time, + surface_list=surface_list, + start_time=start_time, + end_time=end_time, + ) + + # assign the material an rgb value + cell_material_name = cell["material_name"] + if cell_material_name == "water": + rgb = water_rgb + elif cell_material_name == "source": + rgb = source_rgb + else: + rgb = distinct_colors[material_colors_to_generate.index(cell_material_name)] + rgb = [int(rgb[0]), int(rgb[1]), int(rgb[2])] # ngsolve takes rgb as a list + + # if material is missing from the color key, add it + if cell_material_name not in materials_added_to_color_key: + materials_added_to_color_key.append(cell_material_name) + color_key_list.append([rgb, cell_material_name]) + + # add the cell geometry to the visualization + geo.Add(cell_geometry.col(rgb), transparent=True) + + # draw the visualization + geo.Draw() + Redraw() + + return color_key_list + + +# displays the color key to the user +# called by visualize() +def create_color_key(root, color_key_list): + Label(root, text="color key").grid(row=2, column=0) + + # for each material in the color_key_list display + # the material name and corresponding color to the user + for material_index in range(0, len(color_key_list)): + material = color_key_list[material_index] + + # create label for the material name + Label(root, text=str(material[1]) + ":").grid(row=3 + material_index, sticky=W) + + # canvas where color will be displayed + canvas = Canvas(root, width=200, height=len(color_key_list) * 50) + + # switch from rgb to hex for tkinter + rgb = material[0] + rgb = [255 * rgb[0], 255 * rgb[1], 255 * rgb[2]] + colorval = "#{0:02x}{1:02x}{2:02x}".format(rgb[0], rgb[1], rgb[2]) + + # add rectangle to canvas with material color + canvas.create_rectangle(10, 10, 70, 60, fill=colorval) + + # add canvas to the window + canvas.grid(row=3 + material_index, column=1, sticky=E) + + +# triggered when time slider changed +# it redraws the model at the new time +def time_slider_changed(event, start_time, end_time): + draw_Geometry(current_time=float(event), start_time=start_time, end_time=end_time) + + +# creates the time slider +# called by visualize() +def create_time_slider(root, start_time, end_time): + root.title("Time Slider") + + time_label = Label(root, text="Time") + time_label.grid(row=0, column=0, columnspan=2) + + time_scale = Scale( + root, + from_=start_time, + to=end_time, + orient=HORIZONTAL, + tickinterval=1, + command=lambda event: time_slider_changed(event, start_time, end_time), + ) + time_scale.grid(row=1, column=1) + + +# runs the visualization for a model +# start and end times are default zero +# called in input file +def visualize(start_time=0, end_time=0): + import netgen.gui # launches visualiztation window + + color_key_list = draw_Geometry( + current_time=0, start_time=start_time, end_time=end_time + ) + + # Set up tkinter window + root = Tk() + if start_time != end_time: + create_time_slider(root, start_time, end_time) + create_color_key(root, color_key_list) + root.mainloop() # mainloop for tkinter window From 9ad8964dce63a67021961c92117446012bbd4a53 Mon Sep 17 00:00:00 2001 From: "Joanna Piper Morgan (jonsey)" Date: Thu, 21 Dec 2023 17:30:08 -0800 Subject: [PATCH 6/9] altering readme --- README.md | 5 +---- 1 file changed, 1 insertion(+), 4 deletions(-) diff --git a/README.md b/README.md index a814fcde..a86314a5 100644 --- a/README.md +++ b/README.md @@ -10,9 +10,6 @@ MC/DC is a performant, scalable, and machine-portable Python-based Monte Carlo neutron transport software currently developed in the Center for Exascale Monte Carlo Neutron Transport ([CEMeNT](https://cement-psaap.github.io/)). - -This is a test blah blah blah - *Please Note that this project is in the early stages of devlopment (not even an alpha).* That being said enjoy! ## Installation @@ -31,7 +28,7 @@ on the [AZURV1 benchmark](https://inis.iaea.org/search/search.aspx?orig_q=RN:410 ```python import numpy as np -import mcdc test +import mcdc # ============================================================================= # Set model From a69affabeb9422521c3867610836c46e624bbfe1 Mon Sep 17 00:00:00 2001 From: "Joanna Piper Morgan (jonsey)" Date: Thu, 21 Dec 2023 17:38:50 -0800 Subject: [PATCH 7/9] adding new dependencies to install script --- install.sh | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/install.sh b/install.sh index 8134f030..e2c6ae88 100755 --- a/install.sh +++ b/install.sh @@ -6,7 +6,7 @@ pip install -e . # Install MC/DC dependencies, reply "y" to conda prompt -conda install numpy numba matplotlib scipy h5py pytest colorama <<< "y" +conda install numpy numba matplotlib scipy h5py pytest colorama ngsolve distinctipy <<< "y" # Patch numba From c491c068da37fa8cb6cb00d32410220731b15413 Mon Sep 17 00:00:00 2001 From: Sam Pasmann Date: Tue, 2 Jan 2024 23:20:49 -0700 Subject: [PATCH 8/9] fix iQMC MPI bug --- mcdc/main.py | 63 ++++++++++++++++++++------------------------------- mcdc/type_.py | 14 +++++++++--- 2 files changed, 36 insertions(+), 41 deletions(-) diff --git a/mcdc/main.py b/mcdc/main.py index da3e7369..0d0cbb55 100644 --- a/mcdc/main.py +++ b/mcdc/main.py @@ -375,64 +375,51 @@ def prepare(): if input_deck.technique["iQMC"]: # pass in mesh - mcdc["technique"]["iqmc"]["mesh"]["x"] = input_deck.technique["iqmc"]["mesh"][ - "x" - ] - mcdc["technique"]["iqmc"]["mesh"]["y"] = input_deck.technique["iqmc"]["mesh"][ - "y" - ] - mcdc["technique"]["iqmc"]["mesh"]["z"] = input_deck.technique["iqmc"]["mesh"][ - "z" - ] - mcdc["technique"]["iqmc"]["mesh"]["t"] = input_deck.technique["iqmc"]["mesh"][ - "t" - ] + iqmc = mcdc["technique"]["iqmc"] + iqmc["mesh"]["x"] = input_deck.technique["iqmc"]["mesh"]["x"] + iqmc["mesh"]["y"] = input_deck.technique["iqmc"]["mesh"]["y"] + iqmc["mesh"]["z"] = input_deck.technique["iqmc"]["mesh"]["z"] + iqmc["mesh"]["t"] = input_deck.technique["iqmc"]["mesh"]["t"] # pass in score list for name, value in input_deck.technique["iqmc"]["score_list"].items(): - mcdc["technique"]["iqmc"]["score_list"][name] = value + iqmc["score_list"][name] = value # pass in initial tallies for name, value in input_deck.technique["iqmc"]["score"].items(): mcdc["technique"]["iqmc"]["score"][name] = value # LDS generator - mcdc["technique"]["iqmc"]["generator"] = input_deck.technique["iqmc"][ - "generator" - ] + iqmc["generator"] = input_deck.technique["iqmc"]["generator"] # minimum particle weight - mcdc["technique"]["iqmc"]["w_min"] = 1e-16 # / mcdc["setting"]["N_particle"] + iqmc["w_min"] = 1e-13 # variables to generate samples - scramble = mcdc["technique"]["iqmc"]["scramble"] - N_dim = mcdc["technique"]["iqmc"]["N_dim"] - seed = mcdc["technique"]["iqmc"]["seed"] - N = mcdc["setting"]["N_particle"] - size = MPI.COMM_WORLD.Get_size() - rank = MPI.COMM_WORLD.Get_rank() - N_work = math.ceil(N_particle / size) - # how many samples this processor will skip in the LDS - fast_forward = int((rank / size) * N) + scramble = iqmc["scramble"] + N_dim = iqmc["N_dim"] + seed = iqmc["seed"] + + mcdc["mpi_size"] = MPI.COMM_WORLD.Get_size() + mcdc["mpi_rank"] = MPI.COMM_WORLD.Get_rank() + kernel.distribute_work(N_particle, mcdc) + N_work = int(mcdc["mpi_work_size"]) + N_start = int(mcdc["mpi_work_start"]) + # generate LDS if input_deck.technique["iqmc"]["generator"] == "sobol": sampler = qmc.Sobol(d=N_dim, scramble=scramble) # skip the first entry in Sobol sequence because its 0.0 # skip the second because it maps to ux = 0.0 sampler.fast_forward(2) - sampler.fast_forward(fast_forward) - mcdc["technique"]["iqmc"]["lds"] = sampler.random(N_work) + sampler.fast_forward(N_start) + iqmc["lds"] = sampler.random(N_work) if input_deck.technique["iqmc"]["generator"] == "halton": sampler = qmc.Halton(d=N_dim, scramble=scramble, seed=seed) # skip the first entry in Halton sequence because its 0.0 sampler.fast_forward(1) - sampler.fast_forward(fast_forward) - mcdc["technique"]["iqmc"]["lds"] = sampler.random(N_work) + sampler.fast_forward(N_start) + iqmc["lds"] = sampler.random(N_work) if input_deck.technique["iqmc"]["generator"] == "random": - # this chunk of code uses the iqmc_seed to generate a number of - # seeds to be used on each processor - # this way, each processor gets different samples, but if iQMC is run - # several times it will generate the same samples across runs - # 1e6 represents the maximum integer size generated np.random.seed(seed) - seeds = np.random.randint(1e6, size=size) - np.random.seed(seeds[rank]) - mcdc["technique"]["iqmc"]["lds"] = np.random.random((N_work, N_dim)) + seeds = np.random.randint(1e6, size=mcdc["mpi_size"]) + np.random.seed(seeds[mcdc["mpi_rank"]]) + iqmc["lds"] = np.random.random((N_work, N_dim)) # ========================================================================= # Derivative Source Method diff --git a/mcdc/type_.py b/mcdc/type_.py index c0d6e026..76b7409d 100644 --- a/mcdc/type_.py +++ b/mcdc/type_.py @@ -510,8 +510,17 @@ def make_type_technique(N_particle, G, card): iqmc_list += [("mesh", mesh)] # Low-discprenecy sequence - N_work = math.ceil(N_particle / MPI.COMM_WORLD.Get_size()) - iqmc_list += [("lds", float64, (N_work, N_dim))] + size = MPI.COMM_WORLD.Get_size() + rank = MPI.COMM_WORLD.Get_rank() + # Evenly distribute work + work_size = math.floor(N_particle / size) + # Count reminder + rem = N_particle % size + # Assign reminder and update starting index + if rank < rem: + work_size += 1 + + iqmc_list += [("lds", float64, (work_size, N_dim))] iqmc_list += [("fixed_source", float64, (Ng, Nt, Nx, Ny, Nz))] # TODO: make matidx int32 iqmc_list += [("material_idx", int64, (Nt, Nx, Ny, Nz))] @@ -521,7 +530,6 @@ def make_type_technique(N_particle, G, card): iqmc_list += [(("total_source"), float64, (total_size,))] # Scores and shapes - # TODO: add fission power tally scores_shapes = [ ["flux", (Ng, Nt, Nx, Ny, Nz)], ["effective-scattering", (Ng, Nt, Nx, Ny, Nz)], From 8d0ad146dd501b1f551ceca2af29edb5e2016476 Mon Sep 17 00:00:00 2001 From: Sam Pasmann Date: Tue, 2 Jan 2024 23:23:50 -0700 Subject: [PATCH 9/9] add iqmc back to mpi/numba regression test suite --- .github/workflows/tests.yml | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/.github/workflows/tests.yml b/.github/workflows/tests.yml index e224c51f..9032ca6a 100644 --- a/.github/workflows/tests.yml +++ b/.github/workflows/tests.yml @@ -35,8 +35,8 @@ jobs: - name: Regression Test - MPI run: | cd test/regression - python run.py --mpiexec=4 --skip=iqmc* - python run.py --mpiexec=2 --name=iqmc* + python run.py --mpiexec=4 + python run.py --mpiexec=2 - name: Regression Test - Numba run: | cd test/regression @@ -44,7 +44,7 @@ jobs: - name: Regression Test - Numba and MPI run: | cd test/regression - python run.py --mode=numba --mpiexec=4 --skip=iqmc* - python run.py --mode=numba --mpiexec=2 --name=iqmc* + python run.py --mode=numba --mpiexec=4 + python run.py --mode=numba --mpiexec=2 # TODO: Performance test # TODO: GPU test