Merge branch 'main' of https://github.com/ilhamv/MCDC into continuous…

…_energy

.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
@@ -35,16 +35,16 @@ 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
         python run.py --mode=numba
     - 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

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

@@ -216,6 +216,7 @@ def make_card_material(N_nuclide, G=1, J=0):
     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
@@ -245,6 +246,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
 
@@ -257,6 +259,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

@@ -231,6 +231,7 @@ def material(
     chi_d=None,
     speed=None,
     decay=None,
+    name="P",
     sensitivity=False,
     dsm_Np=1.0,
 ):
@@ -344,6 +345,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]
@@ -514,6 +520,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"])
@@ -666,6 +674,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)

mcdc/main.py

@@ -419,64 +419,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_particle = 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

mcdc/type_.py

@@ -738,8 +738,17 @@ def make_type_technique(input_deck):
     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))]
@@ -749,7 +758,6 @@ def make_type_technique(input_deck):
     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)],