Adding visulization (Merging #109)

.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

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()
+"""

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()

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

mcdc/__init__.py

@@ -25,3 +25,4 @@
     reset_cards,
 )
 from mcdc.main import run, prepare
+from mcdc.visualizer import visualize

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])

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)