Merge branch 'dev' into better_tally

.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

docs/source/conf.py

@@ -18,7 +18,7 @@
 # On Read the Docs, need to mock any python packages that would require c
 from unittest.mock import MagicMock
 
-MOCK_MODULES = ["mpi4py", "colorama"]
+MOCK_MODULES = ["mpi4py", "colorama", "mpi4py.util.dtlib"]
 sys.modules.update((mod_name, MagicMock()) for mod_name in MOCK_MODULES)
 
 

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

examples/visulization/godiva-mockup/input.py

@@ -1,5 +1,6 @@
 import numpy as np
 import mcdc, h5py
+from mcdc import visualizer
 
 # =============================================================================
 # Materials
@@ -73,7 +74,7 @@
 
 mcdc.source(x=[-22.0, 22.0], time=[0.0, 5.0], isotropic=True)
 
-mcdc.visualize(
+visualizer.visualize(
     start_time=0, end_time=1, tick_interval=0.1, material_colors={"water": [1, 0, 0]}
 )
 # =============================================================================

mcdc/__init__.py

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

mcdc/adapt.py

@@ -416,6 +416,16 @@ def local_j_array():
     return cuda.local.array(1, type_.j_array)[0]
 
 
+@for_cpu()
+def local_RPN_array():
+    return np.zeros(1, dtype=type_.RPN_array)[0]
+
+
+@for_gpu()
+def local_RPN_array():
+    return cuda.local.array(1, type_.RPN_array)[0]
+
+
 @for_cpu()
 def local_particle():
     return np.zeros(1, dtype=type_.particle)[0]

mcdc/card.py

@@ -1,6 +1,17 @@
 import numpy as np
+import sympy
 
+<<<<<<< HEAD
 from mcdc.constant import INF, SHIFT, PI
+=======
+from mcdc.constant import (
+    BOOL_AND,
+    BOOL_OR,
+    BOOL_NOT,
+    INF,
+    SHIFT,
+)
+>>>>>>> dev
 
 # Get the global variable container
 import mcdc.global_ as global_
@@ -16,7 +27,10 @@ def __str__(self):
         for name in [
             a
             for a in dir(self)
-            if not a.startswith("__") and not callable(getattr(self, a)) and a != "tag"
+            if not a.startswith("__")
+            and not callable(getattr(self, a))
+            and a != "tag"
+            and not a.startswith("_")
         ]:
             text += "  %s : %s\n" % (name, str(getattr(self, name)))
         return text
@@ -56,6 +70,7 @@ def __init__(self, G=1, J=0, name=None):
         self.uq = False
         self.flags = []
         self.distribution = ""
+        self.name = ""
 
 
 class MaterialCard(InputCard):
@@ -93,8 +108,8 @@ def __init__(self, type_):
         # Set card data
         self.ID = None
         self.type = type_
-        self.A = -1
-        self.B = -1
+        self.A = None
+        self.B = None
 
     def __and__(self, other):
         region = RegionCard("intersection")
@@ -122,6 +137,21 @@ def __invert__(self):
         global_.input_deck.regions.append(region)
         return region
 
+    def __str__(self):
+        if self.type == "halfspace":
+            if self.B > 0:
+                return "+s%i" % self.A
+            else:
+                return "-s%i" % self.A
+        elif self.type == "intersection":
+            return "r%i & r%i" % (self.A, self.B)
+        elif self.type == "union":
+            return "r%i | r%i" % (self.A, self.B)
+        elif self.type == "complement":
+            return "~r%i" % (self.A)
+        elif self.type == "all":
+            return "all"
+
 
 class SurfaceCard(InputCard):
     def __init__(self):
@@ -153,23 +183,33 @@ def __init__(self):
         self.N_tally = 0
         self.tally_IDs = []
 
-    def __pos__(self):
+    def _create_halfspace(self, positive):
         region = RegionCard("halfspace")
         region.A = self.ID
-        region.B = 1
+        if positive:
+            region.B = 1
+        else:
+            region.B = -1
+
+        # Check if an identical halfspace region already existed
+        for idx, existing_region in enumerate(global_.input_deck.regions):
+            if (
+                existing_region.type == "halfspace"
+                and region.A == existing_region.A
+                and region.B == existing_region.B
+            ):
+                return global_.input_deck.regions[idx]
+
         # Set ID and push to deck
         region.ID = len(global_.input_deck.regions)
         global_.input_deck.regions.append(region)
         return region
 
+    def __pos__(self):
+        return self._create_halfspace(True)
+
     def __neg__(self):
-        region = RegionCard("halfspace")
-        region.A = self.ID
-        region.B = 0
-        # Set ID and push to deck
-        region.ID = len(global_.input_deck.regions)
-        global_.input_deck.regions.append(region)
-        return region
+        return self._create_halfspace(False)
 
 
 class CellCard(InputCard):
@@ -178,12 +218,87 @@ def __init__(self):
 
         # Set card data
         self.ID = None
-        self.region_ID = -1
+        self.region_ID = None
+        self.region = "all"
         self.fill_type = "material"
-        self.fill_ID = -1
+        self.fill_ID = None
         self.translation = np.array([0.0, 0.0, 0.0])
-        self.N_surface = 0
-        self.surface_ID = np.zeros(0, dtype=int)
+        self.surface_IDs = np.zeros(0, dtype=int)
+        self._region_RPN = []  # Reverse Polish Notation
+
+    def set_region_RPN(self):
+        # Make alias and reset
+        rpn = self._region_RPN
+        rpn.clear()
+
+        # Build RPN based on the assigned region
+        region = global_.input_deck.regions[self.region_ID]
+        stack = [region]
+        while len(stack) > 0:
+            token = stack.pop()
+            if isinstance(token, RegionCard):
+                if token.type == "halfspace":
+                    rpn.append(token.ID)
+                elif token.type == "intersection":
+                    region_A = global_.input_deck.regions[token.A]
+                    region_B = global_.input_deck.regions[token.B]
+                    stack += ["&", region_A, region_B]
+                elif token.type == "union":
+                    region_A = global_.input_deck.regions[token.A]
+                    region_B = global_.input_deck.regions[token.B]
+                    stack += ["|", region_A, region_B]
+                elif token.type == "complement":
+                    region = global_.input_deck.regions[token.A]
+                    stack += ["~", region]
+            else:
+                if token == "&":
+                    rpn.append(BOOL_AND)
+                elif token == "|":
+                    rpn.append(BOOL_OR)
+                elif token == "~":
+                    rpn.append(BOOL_NOT)
+                else:
+                    print_error("Something is wrong with cell RPN creation.")
+
+    def set_region(self):
+        stack = []
+
+        for token in self._region_RPN:
+            if token >= 0:
+                stack.append(token)
+            else:
+                if token == BOOL_AND or token == BOOL_OR:
+                    item_1 = stack.pop()
+                    if isinstance(item_1, int):
+                        item_1 = sympy.symbols(str(global_.input_deck.regions[item_1]))
+
+                    item_2 = stack.pop()
+                    if isinstance(item_2, int):
+                        item_2 = sympy.symbols(str(global_.input_deck.regions[item_2]))
+
+                    if token == BOOL_AND:
+                        stack.append(item_1 & item_2)
+                    else:
+                        stack.append(item_1 | item_2)
+
+                elif token == BOOL_NOT:
+                    item = stack.pop()
+                    if isinstance(item, int):
+                        item = sympy.symbols(str(global_.input_deck.regions[item]))
+                    stack.append(~item)
+
+        self.region = sympy.logic.boolalg.simplify_logic(stack[0])
+
+    def set_surface_IDs(self):
+        surface_IDs = []
+
+        for token in self._region_RPN:
+            if token >= 0:
+                ID = global_.input_deck.regions[token].A
+                if not ID in surface_IDs:
+                    surface_IDs.append(ID)
+
+        self.surface_IDs = np.sort(np.array(surface_IDs))
 
 
 class UniverseCard(InputCard):

mcdc/constant.py

@@ -36,7 +36,12 @@
 REGION_COMPLEMENT = 3
 REGION_ALL = 4
 
-# UNIVERSE
+# Boolean operator
+BOOL_AND = -1
+BOOL_OR = -2
+BOOL_NOT = -3
+
+# Universe
 UNIVERSE_ROOT = 0
 
 # Events