From 85fb12a5a5f8292083fd150e1ad0fecfce99718c Mon Sep 17 00:00:00 2001 From: Ilham Variansyah Date: Thu, 8 Aug 2024 08:03:24 +0700 Subject: [PATCH 1/6] set up the front end --- mcdc/card.py | 129 +++++++++++++++++++++++++++++++++++++++++------ mcdc/constant.py | 7 ++- mcdc/input_.py | 5 ++ 3 files changed, 124 insertions(+), 17 deletions(-) diff --git a/mcdc/card.py b/mcdc/card.py index ecf66167..299fd375 100644 --- a/mcdc/card.py +++ b/mcdc/card.py @@ -1,6 +1,13 @@ import numpy as np +import sympy -from mcdc.constant import INF, SHIFT +from mcdc.constant import ( + BOOL_AND, + BOOL_OR, + BOOL_NOT, + INF, + SHIFT, +) # Get the global variable container import mcdc.global_ as global_ @@ -16,7 +23,7 @@ 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 @@ -88,8 +95,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") @@ -117,6 +124,21 @@ def __invert__(self): global_.input_deck.regions.append(region) return region + def __str__(self): + if self.type == 'halfspace': + if self.B > 0.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): @@ -146,23 +168,32 @@ def __init__(self): self.sensitivity_ID = 0 self.dsm_Np = 1.0 - 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 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): @@ -171,12 +202,78 @@ 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]) + class UniverseCard(InputCard): diff --git a/mcdc/constant.py b/mcdc/constant.py index ab603f52..da60c12a 100644 --- a/mcdc/constant.py +++ b/mcdc/constant.py @@ -20,7 +20,12 @@ REGION_COMPLEMENT = 3 REGION_ALL = 4 -# UNIVERSE +# Boolean operator +BOOL_AND = -1 +BOOL_OR = -2 +BOOL_NOT = -3 + +# Universe UNIVERSE_ROOT = 0 # Events diff --git a/mcdc/input_.py b/mcdc/input_.py index f0a9306f..02079bed 100644 --- a/mcdc/input_.py +++ b/mcdc/input_.py @@ -732,6 +732,11 @@ def cell(region=None, fill=None, translation=(0.0, 0.0, 0.0)): # Assign region card.region_ID = region.ID + # Set region Reverse Polish Notation and region description + if region.type != 'all': + card.set_region_RPN() + card.set_region() + # Assign fill type and ID if fill.tag == "Material": card.fill_type = "material" From 6cba33b4c04f6b1a2c8f69268550708d16e84c91 Mon Sep 17 00:00:00 2001 From: Ilham Variansyah Date: Thu, 8 Aug 2024 08:16:00 +0700 Subject: [PATCH 2/6] add set_surface_IDs to CellCard --- mcdc/card.py | 11 +++++++++++ mcdc/input_.py | 27 +-------------------------- 2 files changed, 12 insertions(+), 26 deletions(-) diff --git a/mcdc/card.py b/mcdc/card.py index 299fd375..797619a5 100644 --- a/mcdc/card.py +++ b/mcdc/card.py @@ -274,6 +274,17 @@ def set_region(self): 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)) + self.N_surface = len(surface_IDs) class UniverseCard(InputCard): diff --git a/mcdc/input_.py b/mcdc/input_.py index 02079bed..848577f2 100644 --- a/mcdc/input_.py +++ b/mcdc/input_.py @@ -750,14 +750,7 @@ def cell(region=None, fill=None, translation=(0.0, 0.0, 0.0)): card.translation[:] = translation # Get all surface IDs - surface_IDs = [] - - region = global_.input_deck.regions[card.region_ID] - get_all_surface_IDs(region, surface_IDs) - surface_IDs = list(set(surface_IDs)) - - card.N_surface = len(surface_IDs) - card.surface_IDs = np.array(surface_IDs) + card.set_surface_IDs() # Add to deck global_.input_deck.cells.append(card) @@ -765,24 +758,6 @@ def cell(region=None, fill=None, translation=(0.0, 0.0, 0.0)): return card -def get_all_surface_IDs(region, surface_IDs): - if region.type == "halfspace": - surface = global_.input_deck.surfaces[region.A] - surface_IDs.append(surface.ID) - - elif region.type in ["intersection", "union"]: - region1 = global_.input_deck.regions[region.A] - region2 = global_.input_deck.regions[region.B] - - get_all_surface_IDs(region1, surface_IDs) - get_all_surface_IDs(region2, surface_IDs) - - elif region.type in ["complement"]: - region1 = global_.input_deck.regions[region.A] - - get_all_surface_IDs(region1, surface_IDs) - - def universe(cells, root=False): """ Define a list of cells as a universe. From bc68a49303a0935fda9fa739ced4b331cf6ccd3e Mon Sep 17 00:00:00 2001 From: Ilham Variansyah Date: Thu, 8 Aug 2024 12:18:52 +0700 Subject: [PATCH 3/6] working on kobayashi-new --- mcdc/adapt.py | 10 +++++++ mcdc/card.py | 4 +-- mcdc/kernel.py | 51 ++++++++++++++++++++++++++++++----- mcdc/main.py | 73 ++++++++++++++++++++++++++++---------------------- mcdc/type_.py | 68 ++++++++++++++++++++-------------------------- 5 files changed, 126 insertions(+), 80 deletions(-) diff --git a/mcdc/adapt.py b/mcdc/adapt.py index 4cb2ac92..d80e6def 100644 --- a/mcdc/adapt.py +++ b/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] diff --git a/mcdc/card.py b/mcdc/card.py index 797619a5..d107cd9f 100644 --- a/mcdc/card.py +++ b/mcdc/card.py @@ -126,7 +126,7 @@ def __invert__(self): def __str__(self): if self.type == 'halfspace': - if self.B > 0.0: + if self.B > 0: return "+s%i"%self.A else: return "-s%i"%self.A @@ -207,7 +207,6 @@ def __init__(self): self.fill_type = "material" self.fill_ID = None self.translation = np.array([0.0, 0.0, 0.0]) - self.N_surface = 0 self.surface_IDs = np.zeros(0, dtype=int) self._region_RPN = [] # Reverse Polish Notation @@ -284,7 +283,6 @@ def set_surface_IDs(self): surface_IDs.append(ID) self.surface_IDs = np.sort(np.array(surface_IDs)) - self.N_surface = len(surface_IDs) class UniverseCard(InputCard): diff --git a/mcdc/kernel.py b/mcdc/kernel.py index 4e689ee7..72461eb3 100644 --- a/mcdc/kernel.py +++ b/mcdc/kernel.py @@ -1708,8 +1708,40 @@ def split_particle(P): @njit def cell_check(P, cell, trans, mcdc): - region = mcdc["regions"][cell["region_ID"]] - return region_check(P, region, trans, mcdc) + # Access RPN data + idx = cell['region_data_idx'] + N_token = mcdc['cell_region_data'][idx] + + # Create local value array + value_struct = adapt.local_RPN_array() + value = value_struct['values'] + N_value = 0 + + # March forward through RPN tokens + idx += 1 + idx_end = idx + N_token + while idx < idx_end: + token = mcdc['cell_region_data'][idx] + + if token >= 0: + surface = mcdc['surfaces'][token] + value[N_value] = surface_evaluate(P, surface, trans) > 0.0 + N_value += 1 + + elif token == BOOL_NOT: + value[N_value - 1] = not value[N_value - 1] + + elif token == BOOL_AND: + value[N_value - 2] = value[N_value - 2] & value[N_value - 1] + N_value -= 1 + + elif token == BOOL_OR: + value[N_value - 2] = value[N_value - 2] | value[N_value - 1] + N_value -= 1 + + idx += 1 + + return value[0] @njit @@ -1864,7 +1896,6 @@ def surface_shift(P, surface, trans, mcdc): dot = ux * nx + uy * ny + uz * nz + J1 / v else: dot = ux * nx + uy * ny + uz * nz - if dot > 0.0: P["x"] += shift_x P["y"] += shift_y @@ -2531,7 +2562,6 @@ def move_to_event(P, mcdc): # Also set particle material and speed d_boundary, event = distance_to_boundary(P, mcdc) - # print(P["material_ID"]) # Distance to tally mesh d_mesh = INF if mcdc["cycle_active"]: @@ -2693,13 +2723,22 @@ def distance_to_nearest_surface(P, cell, trans, mcdc): surface_ID = -1 surface_move = False - for i in range(cell["N_surface"]): - surface = mcdc["surfaces"][cell["surface_IDs"][i]] + # Access cell surface data + idx = cell['surface_data_idx'] + N_surface = mcdc['cell_surface_data'][idx] + + # Iterate over all surfaces + idx += 1 + idx_end = idx + N_surface + while idx < idx_end: + candidate_surface_ID = mcdc['cell_surface_data'][idx] + surface = mcdc["surfaces"][candidate_surface_ID] d, sm = surface_distance(P, surface, trans, mcdc) if d < distance: distance = d surface_ID = surface["ID"] surface_move = sm + idx += 1 return distance, surface_ID, surface_move diff --git a/mcdc/main.py b/mcdc/main.py index ffa9d264..2803540f 100644 --- a/mcdc/main.py +++ b/mcdc/main.py @@ -337,6 +337,24 @@ def prepare(): root_universe.cell_IDs[i] = cell.ID input_deck.universes[0] = root_universe + # ========================================================================= + # Prepare cell region RPN (Reverse Polish Notation) + # - Replace halfspace region ID with its surface and insert + # complement operator if the sense is negative. + # ========================================================================= + + for cell in input_deck.cells: + i = 0 + while i < len(cell._region_RPN): + token = cell._region_RPN[i] + if token >= 0: + surface_ID = input_deck.regions[token].A + sense = input_deck.regions[token].B + cell._region_RPN[i] = surface_ID + if sense < 0: + cell._region_RPN.insert(i + 1, BOOL_NOT) + i += 1 + # ========================================================================= # Adapt kernels # ========================================================================= @@ -352,8 +370,6 @@ def prepare(): type_.make_type_nuclide(input_deck) type_.make_type_material(input_deck) type_.make_type_surface(input_deck) - type_.make_type_region() - type_.make_type_cell(input_deck) type_.make_type_universe(input_deck) type_.make_type_lattice(input_deck) type_.make_type_source(input_deck) @@ -370,6 +386,7 @@ def prepare(): type_.make_type_translate(input_deck) type_.make_type_group_array(input_deck) type_.make_type_j_array(input_deck) + type_.make_type_RPN_array(input_deck) # ========================================================================= # Create the global variable container @@ -519,37 +536,16 @@ def prepare(): N = len(getattr(input_deck.surfaces[i], name)) mcdc["surfaces"][i][name][:N] = getattr(input_deck.surfaces[i], name) - # ========================================================================= - # Regions - # ========================================================================= - - N_region = len(input_deck.regions) - for i in range(N_region): - for name in type_.region.names: - if name not in ["type"]: - copy_field(mcdc["regions"][i], input_deck.regions[i], name) - - # Type - if input_deck.regions[i].type == "halfspace": - mcdc["regions"][i]["type"] = REGION_HALFSPACE - elif input_deck.regions[i].type == "intersection": - mcdc["regions"][i]["type"] = REGION_INTERSECTION - elif input_deck.regions[i].type == "union": - mcdc["regions"][i]["type"] = REGION_UNION - elif input_deck.regions[i].type == "complement": - mcdc["regions"][i]["type"] = REGION_COMPLEMENT - elif input_deck.regions[i].type == "all": - mcdc["regions"][i]["type"] = REGION_ALL - # ========================================================================= # Cells # ========================================================================= N_cell = len(input_deck.cells) + surface_data_idx = 0 + region_data_idx = 0 for i in range(N_cell): - for name in type_.cell.names: - if name not in ["fill_type", "surface_IDs"]: - copy_field(mcdc["cells"][i], input_deck.cells[i], name) + for name in ['ID', 'fill_ID', 'translation']: + copy_field(mcdc["cells"][i], input_deck.cells[i], name) # Fill type if input_deck.cells[i].fill_type == "material": @@ -559,10 +555,19 @@ def prepare(): elif input_deck.cells[i].fill_type == "lattice": mcdc["cells"][i]["fill_type"] = FILL_LATTICE - # Variables with possible different sizes - for name in ["surface_IDs"]: - N = mcdc["cells"][i]["N_surface"] - mcdc["cells"][i][name][:N] = getattr(input_deck.cells[i], name) + # Surface data + mcdc["cells"][i]['surface_data_idx'] = surface_data_idx + N_surface = len(input_deck.cells[i].surface_IDs) + mcdc['cell_surface_data'][surface_data_idx] = N_surface + mcdc['cell_surface_data'][surface_data_idx + 1: surface_data_idx + N_surface + 1] = input_deck.cells[i].surface_IDs + surface_data_idx += N_surface + 1 + + # Region data + mcdc["cells"][i]['region_data_idx'] = region_data_idx + N_RPN = len(input_deck.cells[i]._region_RPN) + mcdc['cell_region_data'][region_data_idx] = N_RPN + mcdc['cell_region_data'][region_data_idx + 1: region_data_idx + N_RPN + 1] = input_deck.cells[i]._region_RPN + region_data_idx += N_RPN + 1 # ========================================================================= # Universes @@ -981,7 +986,11 @@ def card_to_h5group(card, group): elif value is None: next else: - group[name] = value + if name not in ['region']: + group[name] = value + + elif name == 'region': + group[name] = str(value) def dictlist_to_h5group(dictlist, input_group, name): diff --git a/mcdc/type_.py b/mcdc/type_.py index 3dfb9b2c..b17c38d6 100644 --- a/mcdc/type_.py +++ b/mcdc/type_.py @@ -34,8 +34,6 @@ nuclide = None material = None surface = None -region = None -cell = None universe = None lattice = None source = None @@ -47,6 +45,8 @@ translate = None group_array = None j_array = None +RPN_array = None + global_ = None @@ -513,46 +513,21 @@ def make_type_surface(input_deck): ) -# ============================================================================== -# Region -# ============================================================================== - - -def make_type_region(): - global region - - region = into_dtype( - [ - ("ID", int64), - ("type", int64), - ("A", int64), - ("B", int64), - ] - ) - - # ============================================================================== # Cell # ============================================================================== -def make_type_cell(input_deck): - global cell - - # Maximum number of surfaces per cell - Nmax_surface = max([cell.N_surface for cell in input_deck.cells]) - - cell = into_dtype( - [ - ("ID", int64), - ("region_ID", int64), - ("fill_type", int64), - ("fill_ID", int64), - ("translation", float64, (3,)), - ("N_surface", int64), - ("surface_IDs", int64, (Nmax_surface,)), - ] - ) +cell = into_dtype( + [ + ("ID", int64), + ("fill_type", int64), + ("fill_ID", int64), + ("translation", float64, (3,)), + ("surface_data_idx", int64), + ("region_data_idx", int64), + ] +) # ============================================================================== @@ -1233,12 +1208,18 @@ def make_type_global(input_deck): N_nuclide = len(input_deck.nuclides) N_material = len(input_deck.materials) N_surface = len(input_deck.surfaces) - N_region = len(input_deck.regions) N_cell = len(input_deck.cells) N_source = len(input_deck.sources) N_universe = len(input_deck.universes) N_lattice = len(input_deck.lattices) + # Cell data sizes + N_cell_surface = 0 + N_cell_region = 0 + for cell_ in input_deck.cells: + N_cell_surface += 1 + len(cell_.surface_IDs) + N_cell_region += 1 + len(cell_._region_RPN) + # Simulation parameters N_particle = input_deck.setting["N_particle"] N_precursor = input_deck.setting["N_precursor"] @@ -1298,8 +1279,9 @@ def make_type_global(input_deck): ("nuclides", nuclide, (N_nuclide,)), ("materials", material, (N_material,)), ("surfaces", surface, (N_surface,)), - ("regions", region, (N_region,)), ("cells", cell, (N_cell,)), + ("cell_surface_data", int64, (N_cell_surface,)), + ("cell_region_data", int64, (N_cell_region,)), ("universes", universe, (N_universe,)), ("lattices", lattice, (N_lattice,)), ("sources", source, (N_source,)), @@ -1385,6 +1367,14 @@ def make_type_j_array(input_deck): j_array = into_dtype([("values", float64, (J,))]) +def make_type_RPN_array(input_deck): + global RPN_array + N_max = 0 + for cell_ in input_deck.cells: + N = np.sum(np.array(cell_._region_RPN) >= 0.0) + N_max = max(N_max, N) + RPN_array = into_dtype([("values", bool_, (N_max,))]) + def make_type_mesh(card): Nx = len(card["x"]) - 1 Ny = len(card["y"]) - 1 From b01a039759b6ef2d7554fb3be81b72a15ce741c7 Mon Sep 17 00:00:00 2001 From: Ilham Variansyah Date: Thu, 8 Aug 2024 15:48:32 +0700 Subject: [PATCH 4/6] debug halfspace checking --- mcdc/card.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/mcdc/card.py b/mcdc/card.py index d107cd9f..a5670127 100644 --- a/mcdc/card.py +++ b/mcdc/card.py @@ -179,7 +179,7 @@ def _create_halfspace(self, positive): # Check if an identical halfspace region already existed for idx, existing_region in enumerate(global_.input_deck.regions): - if region.A == existing_region.A and region.B == existing_region.B: + 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 From bc81119611f867df38bd198954f9e26e020f53ff Mon Sep 17 00:00:00 2001 From: Ilham Variansyah Date: Thu, 8 Aug 2024 15:51:12 +0700 Subject: [PATCH 5/6] add sympy to pyproject.toml --- pyproject.toml | 1 + 1 file changed, 1 insertion(+) diff --git a/pyproject.toml b/pyproject.toml index 3dc378b2..4fed70cd 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -56,6 +56,7 @@ dependencies = [ "h5py", "colorama", "black", + "sympy", ] [project.optional-dependencies] docs = ["sphinx==7.2.6", "furo", "sphinx_toolbox"] From f793529b6b391f10d5e0a86cf4563d3e4b7be50a Mon Sep 17 00:00:00 2001 From: Ilham Variansyah Date: Thu, 8 Aug 2024 15:53:40 +0700 Subject: [PATCH 6/6] back in black --- mcdc/card.py | 58 +++++++++++++++++++++++++++----------------------- mcdc/input_.py | 2 +- mcdc/kernel.py | 16 +++++++------- mcdc/main.py | 22 +++++++++++-------- mcdc/type_.py | 1 + 5 files changed, 54 insertions(+), 45 deletions(-) diff --git a/mcdc/card.py b/mcdc/card.py index a5670127..fa53a1fa 100644 --- a/mcdc/card.py +++ b/mcdc/card.py @@ -23,7 +23,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" and not a.startswith("_") + 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 @@ -125,18 +128,18 @@ def __invert__(self): return region def __str__(self): - if self.type == 'halfspace': + if self.type == "halfspace": if self.B > 0: - return "+s%i"%self.A + 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 "-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" @@ -168,7 +171,6 @@ def __init__(self): self.sensitivity_ID = 0 self.dsm_Np = 1.0 - def _create_halfspace(self, positive): region = RegionCard("halfspace") region.A = self.ID @@ -179,7 +181,11 @@ def _create_halfspace(self, positive): # 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: + 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 @@ -187,11 +193,9 @@ def _create_halfspace(self, positive): global_.input_deck.regions.append(region) return region - def __pos__(self): return self._create_halfspace(True) - def __neg__(self): return self._create_halfspace(False) @@ -208,7 +212,7 @@ def __init__(self): self.fill_ID = None self.translation = np.array([0.0, 0.0, 0.0]) self.surface_IDs = np.zeros(0, dtype=int) - self._region_RPN = [] # Reverse Polish Notation + self._region_RPN = [] # Reverse Polish Notation def set_region_RPN(self): # Make alias and reset @@ -221,25 +225,25 @@ def set_region_RPN(self): while len(stack) > 0: token = stack.pop() if isinstance(token, RegionCard): - if token.type == 'halfspace': + if token.type == "halfspace": rpn.append(token.ID) - elif token.type == 'intersection': + 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': + 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': + stack += ["|", region_A, region_B] + elif token.type == "complement": region = global_.input_deck.regions[token.A] - stack += ['~', region] + stack += ["~", region] else: - if token == '&': + if token == "&": rpn.append(BOOL_AND) - elif token == '|': + elif token == "|": rpn.append(BOOL_OR) - elif token == '~': + elif token == "~": rpn.append(BOOL_NOT) else: print_error("Something is wrong with cell RPN creation.") @@ -277,7 +281,7 @@ def set_surface_IDs(self): surface_IDs = [] for token in self._region_RPN: - if token >=0: + if token >= 0: ID = global_.input_deck.regions[token].A if not ID in surface_IDs: surface_IDs.append(ID) diff --git a/mcdc/input_.py b/mcdc/input_.py index 848577f2..02c14c77 100644 --- a/mcdc/input_.py +++ b/mcdc/input_.py @@ -733,7 +733,7 @@ def cell(region=None, fill=None, translation=(0.0, 0.0, 0.0)): card.region_ID = region.ID # Set region Reverse Polish Notation and region description - if region.type != 'all': + if region.type != "all": card.set_region_RPN() card.set_region() diff --git a/mcdc/kernel.py b/mcdc/kernel.py index 72461eb3..a02f321f 100644 --- a/mcdc/kernel.py +++ b/mcdc/kernel.py @@ -1709,22 +1709,22 @@ def split_particle(P): @njit def cell_check(P, cell, trans, mcdc): # Access RPN data - idx = cell['region_data_idx'] - N_token = mcdc['cell_region_data'][idx] + idx = cell["region_data_idx"] + N_token = mcdc["cell_region_data"][idx] # Create local value array value_struct = adapt.local_RPN_array() - value = value_struct['values'] + value = value_struct["values"] N_value = 0 # March forward through RPN tokens idx += 1 idx_end = idx + N_token while idx < idx_end: - token = mcdc['cell_region_data'][idx] + token = mcdc["cell_region_data"][idx] if token >= 0: - surface = mcdc['surfaces'][token] + surface = mcdc["surfaces"][token] value[N_value] = surface_evaluate(P, surface, trans) > 0.0 N_value += 1 @@ -2724,14 +2724,14 @@ def distance_to_nearest_surface(P, cell, trans, mcdc): surface_move = False # Access cell surface data - idx = cell['surface_data_idx'] - N_surface = mcdc['cell_surface_data'][idx] + idx = cell["surface_data_idx"] + N_surface = mcdc["cell_surface_data"][idx] # Iterate over all surfaces idx += 1 idx_end = idx + N_surface while idx < idx_end: - candidate_surface_ID = mcdc['cell_surface_data'][idx] + candidate_surface_ID = mcdc["cell_surface_data"][idx] surface = mcdc["surfaces"][candidate_surface_ID] d, sm = surface_distance(P, surface, trans, mcdc) if d < distance: diff --git a/mcdc/main.py b/mcdc/main.py index 2803540f..d47f7e91 100644 --- a/mcdc/main.py +++ b/mcdc/main.py @@ -544,7 +544,7 @@ def prepare(): surface_data_idx = 0 region_data_idx = 0 for i in range(N_cell): - for name in ['ID', 'fill_ID', 'translation']: + for name in ["ID", "fill_ID", "translation"]: copy_field(mcdc["cells"][i], input_deck.cells[i], name) # Fill type @@ -556,17 +556,21 @@ def prepare(): mcdc["cells"][i]["fill_type"] = FILL_LATTICE # Surface data - mcdc["cells"][i]['surface_data_idx'] = surface_data_idx + mcdc["cells"][i]["surface_data_idx"] = surface_data_idx N_surface = len(input_deck.cells[i].surface_IDs) - mcdc['cell_surface_data'][surface_data_idx] = N_surface - mcdc['cell_surface_data'][surface_data_idx + 1: surface_data_idx + N_surface + 1] = input_deck.cells[i].surface_IDs + mcdc["cell_surface_data"][surface_data_idx] = N_surface + mcdc["cell_surface_data"][ + surface_data_idx + 1 : surface_data_idx + N_surface + 1 + ] = input_deck.cells[i].surface_IDs surface_data_idx += N_surface + 1 # Region data - mcdc["cells"][i]['region_data_idx'] = region_data_idx + mcdc["cells"][i]["region_data_idx"] = region_data_idx N_RPN = len(input_deck.cells[i]._region_RPN) - mcdc['cell_region_data'][region_data_idx] = N_RPN - mcdc['cell_region_data'][region_data_idx + 1: region_data_idx + N_RPN + 1] = input_deck.cells[i]._region_RPN + mcdc["cell_region_data"][region_data_idx] = N_RPN + mcdc["cell_region_data"][region_data_idx + 1 : region_data_idx + N_RPN + 1] = ( + input_deck.cells[i]._region_RPN + ) region_data_idx += N_RPN + 1 # ========================================================================= @@ -986,10 +990,10 @@ def card_to_h5group(card, group): elif value is None: next else: - if name not in ['region']: + if name not in ["region"]: group[name] = value - elif name == 'region': + elif name == "region": group[name] = str(value) diff --git a/mcdc/type_.py b/mcdc/type_.py index b17c38d6..43b0e333 100644 --- a/mcdc/type_.py +++ b/mcdc/type_.py @@ -1375,6 +1375,7 @@ def make_type_RPN_array(input_deck): N_max = max(N_max, N) RPN_array = into_dtype([("values", bool_, (N_max,))]) + def make_type_mesh(card): Nx = len(card["x"]) - 1 Ny = len(card["y"]) - 1