Merge pull request #215 from ilhamv/bug_fix

Bug fix, continuous energy test, and file-based source test

.github/workflows/numba_gpu.yml

@@ -34,4 +34,4 @@ jobs:
     - name: Regression Test - Numba
       run: |
           cd test/regression
-          python run.py --mode=numba --skip=iqmc*
+          python run.py --mode=numba

.github/workflows/verification_man_mpi_numba.yml

@@ -32,4 +32,4 @@ jobs:
     - name: Verification Tests - Numba and MPI
       run: |
           cd test/verification/analytic
-          python run.py --mode=numba --mpiexec=2 --skip=iqmc*
+          python run.py --mode=numba --mpiexec=2

.gitignore

@@ -55,3 +55,7 @@ docs/build
 docs/source/pythonapi/generated/
 
 *.csv
+
+# Regression tests
+dummy_nuclide.h5
+source_particles.h5

examples/fixed_source/slab_ce/input.py

@@ -61,9 +61,9 @@
 s4 = mcdc.surface("plane-x", x=2.0, bc="reflective")
 
 # Set cells
-mcdc.cell([+s1, -s2], uo2)
-mcdc.cell([+s2, -s3], h2o)
-mcdc.cell([+s3, -s4], b4c)
+mcdc.cell(+s1 & -s2, uo2)
+mcdc.cell(+s2 & -s3, h2o)
+mcdc.cell(+s3 & -s4, b4c)
 
 # =============================================================================
 # Set source

mcdc/__init__.py

@@ -22,6 +22,8 @@
     uq,
     reset,
     domain_decomposition,
+    make_particle_bank,
+    save_particle_bank,
 )
 from mcdc.main import run, prepare
 from mcdc.visualizer import visualize

mcdc/card.py

@@ -50,6 +50,7 @@ def __init__(self, G=1, J=0):
         self.uq = False
         self.flags = []
         self.distribution = ""
+        self.name = ""
 
 
 class MaterialCard(InputCard):

mcdc/input_.py

@@ -322,26 +322,51 @@ def material(
         # Set ID
         card.ID = len(global_.input_deck.materials)
 
+        # Default values
+        card.J = 6
+
         # Set the nuclides
         for i in range(N_nuclide):
             nuc_name = nuclides[i][0]
             density = nuclides[i][1]
+
+            # Create nuclide card if not defined yet
             if not nuclide_registered(nuc_name):
                 nuc_card = NuclideCard()
+                nuc_card.name = nuc_name
+
+                # Set ID
+                nuc_card.ID = len(global_.input_deck.nuclides)
 
+                # Default values
+                nuc_card.J = 6
+
+                # Check if the nuclide is available in the nuclear data library
                 dir_name = os.getenv("MCDC_XSLIB")
                 if dir_name == None:
                     print_error(
-                        "Continuous energy data directory not configured \n       see https://cement-psaapgithubio.readthedocs.io/en/latest/install.html#configuring-continuous-energy-library \n"
+                        "Continuous energy data directory not configured \n       "
+                        "see https://cement-psaapgithubio.readthedocs.io/en/latest"
+                        "/install.html#configuring-continuous-energy-library \n"
                     )
+
+                # Fissionable flag
                 with h5py.File(dir_name + "/" + nuc_name + ".h5", "r") as f:
                     if max(f["fission"][:]) > 0.0:
                         nuc_card.fissionable = True
+                        card.fissionable = True
+
+                # Add to deck
+                global_.input_deck.nuclides.append(nuc_card)
             else:
                 nuc_card = get_nuclide(nuc_name)
+
             card.nuclide_IDs[i] = nuc_card.ID
             card.nuclide_densities[i] = density
 
+        # Add to deck
+        global_.input_deck.materials.append(card)
+
         return card
 
     # Nuclide and group sizes
@@ -1888,14 +1913,14 @@ def append_card(delta_card, global_tag):
 
 def nuclide_registered(name):
     for card in global_.input_deck.nuclides:
-        if name == card["name"]:
+        if name == card.name:
             return True
     return False
 
 
 def get_nuclide(name):
     for card in global_.input_deck.nuclides:
-        if name == card["name"]:
+        if name == card.name:
             return card
 
 
@@ -1933,6 +1958,41 @@ def check_requirement(label, kw, required):
         print_error("Parameters " + missing + " are required for" + label)
 
 
+def make_particle_bank(size):
+    struct = [
+        ("x", np.float64),
+        ("y", np.float64),
+        ("z", np.float64),
+        ("t", np.float64),
+        ("ux", np.float64),
+        ("uy", np.float64),
+        ("uz", np.float64),
+        ("g", np.uint64),
+        ("E", np.float64),
+        ("w", np.float64),
+        ("sensitivity_ID", np.int64),
+        ("rng_seed", np.uint64),
+    ]
+    iqmc_struct = [("w", np.float64, (1,))]
+    struct += [("iqmc", iqmc_struct)]
+
+    bank = np.zeros(size, dtype=np.dtype(struct))
+
+    # Set default values
+    for i in range(size):
+        bank[i]["ux"] = 1.0
+        bank[i]["w"] = 1.0
+        bank[i]["rng_seed"] = 1
+
+    return bank
+
+
+def save_particle_bank(bank, name):
+    with h5py.File(name + ".h5", "w") as f:
+        f.create_dataset("particles", data=bank[:])
+        f.create_dataset("particles_size", data=len(bank[:]))
+
+
 # ==============================================================================
 # Reset
 # ==============================================================================

mcdc/iqmc/iqmc_loop.py

@@ -79,8 +79,8 @@ def iqmc_step_particle(P, prog):
         kernel.surface_crossing(P, prog)
         if event & EVENT_DOMAIN:
             if not (
-                mcdc["surfaces"][P["surface_ID"]]["reflective"]
-                or mcdc["surfaces"][P["surface_ID"]]["vacuum"]
+                mcdc["surfaces"][P["surface_ID"]]["BC"] == BC_REFLECTIVE
+                or mcdc["surfaces"][P["surface_ID"]]["BC"] == BC_VACUUM
             ):
                 kernel.domain_crossing(P, mcdc)
 

mcdc/main.py

@@ -438,7 +438,7 @@ def prepare():
         # CE data (load data from XS library)
         dir_name = os.getenv("MCDC_XSLIB")
         if mode_CE:
-            nuc_name = input_deck.nuclides[i]["name"]
+            nuc_name = input_deck.nuclides[i].name
             with h5py.File(dir_name + "/" + nuc_name + ".h5", "r") as f:
                 # Atomic weight ratio
                 mcdc["nuclides"][i]["A"] = f["A"][()]
@@ -885,22 +885,29 @@ def prepare():
 
     # =========================================================================
     # Source file
+    #   TODO: Use parallel h5py
     # =========================================================================
 
-    if mcdc["setting"]["source_file"]:
-        with h5py.File(mcdc["setting"]["source_file_name"], "r") as f:
-            # Get source particle size
-            N_particle = f["particles_size"][()]
-
-            # Redistribute work
-            kernel.distribute_work(N_particle, mcdc)
-            N_local = mcdc["mpi_work_size"]
-            start = mcdc["mpi_work_start"]
-            end = start + N_local
-
-            # Add particles to source bank
-            mcdc["bank_source"]["particles"][:N_local] = f["particles"][start:end]
-            mcdc["bank_source"]["size"] = N_local
+    # All ranks, take turn
+    for i in range(mcdc["mpi_size"]):
+        if mcdc["mpi_rank"] == i:
+            if mcdc["setting"]["source_file"]:
+                with h5py.File(mcdc["setting"]["source_file_name"], "r") as f:
+                    # Get source particle size
+                    N_particle = f["particles_size"][()]
+
+                    # Redistribute work
+                    kernel.distribute_work(N_particle, mcdc)
+                    N_local = mcdc["mpi_work_size"]
+                    start = mcdc["mpi_work_start"]
+                    end = start + N_local
+
+                    # Add particles to source bank
+                    mcdc["bank_source"]["particles"][:N_local] = f["particles"][
+                        start:end
+                    ]
+                    mcdc["bank_source"]["size"] = N_local
+        MPI.COMM_WORLD.Barrier()
 
     # =========================================================================
     # IC file

mcdc/type_.py

@@ -315,7 +315,7 @@ def make_type_nuclide(input_deck):
 
         dir_name = os.getenv("MCDC_XSLIB")
         for nuc in input_deck.nuclides:
-            with h5py.File(dir_name + "/" + nuc["name"] + ".h5", "r") as f:
+            with h5py.File(dir_name + "/" + nuc.name + ".h5", "r") as f:
                 NE_xs = max(NE_xs, len(f["E_xs"][:]))
                 NE_nu_p = max(NE_nu_p, len(f["E_nu_p"][:]))
                 NE_nu_d = max(NE_nu_d, len(f["E_nu_d"][:]))
@@ -624,7 +624,7 @@ def make_type_source(input_deck):
     if mode_CE:
         G = 1
         # Maximum number of data point in energy pdf
-        Nmax_E = max([source["energy"].shape[1] for source in input_deck.sources])
+        Nmax_E = max([source.energy.shape[1] for source in input_deck.sources])
     if mode_MG:
         G = input_deck.materials[0].G
         Nmax_E = 2