diff --git a/openmc/cell.py b/openmc/cell.py index 06a80e173c..e01cb86137 100644 --- a/openmc/cell.py +++ b/openmc/cell.py @@ -422,12 +422,18 @@ class Cell(IDManagerMixin): return universes - def clone(self, memo=None): + def clone(self, clone_materials=True, clone_regions=True, memo=None): """Create a copy of this cell with a new unique ID, and clones the cell's region and fill. Parameters ---------- + clone_materials : bool + Whether to create separate copies of the materials filling cells + contained in this cell, or the material filling this cell. + clone_regions : bool + Whether to create separate copies of the regions bounding cells + contained in this cell, and the region bounding this cell. memo : dict or None A nested dictionary of previously cloned objects. This parameter is used internally and should not be specified by the user. @@ -456,13 +462,25 @@ class Cell(IDManagerMixin): self._paths = paths if self.region is not None: - clone.region = self.region.clone(memo) + if clone_regions: + clone.region = self.region.clone(memo) + else: + clone.region = self.region if self.fill is not None: if self.fill_type == 'distribmat': - clone.fill = [fill.clone(memo) if fill is not None else None - for fill in self.fill] + if not clone_materials: + clone.fill = self.fill + else: + clone.fill = [fill.clone(memo) if fill is not None else + None for fill in self.fill] + elif self.fill_type == 'material': + if not clone_materials: + clone.fill = self.fill + else: + clone.fill = self.fill.clone(memo) else: - clone.fill = self.fill.clone(memo) + clone.fill = self.fill.clone(clone_materials, + clone_regions, memo) # Memoize the clone memo[self] = clone diff --git a/openmc/geometry.py b/openmc/geometry.py index b911b01654..205fe8face 100644 --- a/openmc/geometry.py +++ b/openmc/geometry.py @@ -386,7 +386,8 @@ class Geometry(object): surfaces = OrderedDict() for cell in self.get_all_cells().values(): - surfaces = cell.region.get_surfaces(surfaces) + if cell.region is not None: + surfaces = cell.region.get_surfaces(surfaces) return surfaces def get_redundant_surfaces(self): diff --git a/openmc/lattice.py b/openmc/lattice.py index 71176c8ba8..c68637fa34 100644 --- a/openmc/lattice.py +++ b/openmc/lattice.py @@ -7,6 +7,8 @@ from numbers import Real from xml.etree import ElementTree as ET import numpy as np +import warnings +import types import openmc.checkvalue as cv import openmc @@ -454,12 +456,18 @@ class Lattice(IDManagerMixin, metaclass=ABCMeta): return [] return [(self, idx)] + u.find(p) - def clone(self, memo=None): + def clone(self, clone_materials=True, clone_regions=True, memo=None): """Create a copy of this lattice with a new unique ID, and clones all universes within this lattice. Parameters ---------- + clone_materials : bool + Whether to create separate copies of the materials filling cells + contained in this lattice and its outer universe. + clone_regions : bool + Whether to create separate copies of the regions bounding cells + contained in this lattice and its outer universe. memo : dict or None A nested dictionary of previously cloned objects. This parameter is used internally and should not be specified by the user. @@ -480,19 +488,23 @@ class Lattice(IDManagerMixin, metaclass=ABCMeta): clone.id = None if self.outer is not None: - clone.outer = self.outer.clone(memo) + clone.outer = self.outer.clone(clone_materials, clone_regions, + memo) # Assign universe clones to the lattice clone for i in self.indices: if isinstance(self, RectLattice): - clone.universes[i] = self.universes[i].clone(memo) + clone.universes[i] = self.universes[i].clone( + clone_materials, clone_regions, memo) else: if self.ndim == 2: clone.universes[i[0]][i[1]] = \ - self.universes[i[0]][i[1]].clone(memo) + self.universes[i[0]][i[1]].clone(clone_materials, + clone_regions, memo) else: clone.universes[i[0]][i[1]][i[2]] = \ - self.universes[i[0]][i[1]][i[2]].clone(memo) + self.universes[i[0]][i[1]][i[2]].clone( + clone_materials, clone_regions, memo) # Memoize the clone memo[self] = clone @@ -759,6 +771,245 @@ class RectLattice(Lattice): 0 <= idx[1] < self.shape[1] and 0 <= idx[2] < self.shape[2]) + def discretize(self, strategy="degenerate", + universes_to_ignore=[], + materials_to_clone=[], + lattice_neighbors=[], key=lambda univ: univ.id): + """Discretize the lattice with either a degenerate or a local neighbor + symmetry strategy + + 'Degenerate' clones every universe in the lattice, thus making them all + uniquely defined. This is typically required if depletion or thermal + hydraulics will make every universe's environment unique. + + 'Local neighbor symmetry' groups universes with similar neighborhoods. + These clusters of cells and materials provide increased convergence + speed to multi-group cross sections tallies. The user can specify + the lattice's neighbors to discriminate between two sides of a + lattice for example. + + Parameters + ---------- + strategy : {'degenerate', 'lns'} + Which strategy to adopt when discretizing the lattice + universes_to_ignore : Iterable of Universe + Lattice universes that need not be discretized + materials_to_clone : Iterable of Material + List of materials that should be cloned when discretizing + lattice_neighbors : Iterable of Universe + List of the lattice's neighbors. By default, if present, the + lattice outer universe will be used. The neighbors should be + ordered as follows [top left, top, top right, left, right, + bottom left, bottom, bottom right] + key : function + Function of argument a universe that is used to extract a + comparison key. This function will be called on each universe's + neighbors in the lattice to form a neighbor pattern. This pattern + is then used to identify unique neighbor symmetries. + """ + + # Check routine inputs + if self.ndim != 2: + raise NotImplementedError("LNS discretization is not implemented " + "for 1D and 3D lattices") + + cv.check_value('strategy', strategy, ('degenerate', 'lns')) + cv.check_type('universes_to_ignore', universes_to_ignore, Iterable, + openmc.Universe) + cv.check_type('materials_to_clone', materials_to_clone, Iterable, + openmc.Material) + cv.check_type('lattice_neighbors', lattice_neighbors, Iterable, + openmc.Universe) + cv.check_value('number of lattice_neighbors', len(lattice_neighbors), + (0, 8)) + cv.check_type('key', key, types.FunctionType) + + # Use outer universe if neighbors are missing and outer is defined + if self.outer is not None and len(lattice_neighbors) == 0: + lattice_neighbors = [key(self.outer) for i in range(8)] + elif len(lattice_neighbors) == 8: + lattice_neighbors = [key(universe) for universe in + lattice_neighbors] + + # Dictionary that will keep track of where each pattern appears, how + # it was rotated and/or symmetrized + patterns = {} + + # Initialize pattern array + pattern = np.empty(shape=(3, 3), dtype=type(key(self.universes[0][0]))) + + # Define an auxiliary function that returns a universe's neighbors + # that are outside the lattice + def find_edge_neighbors(pattern, i, j): + + # If no neighbors have been specified, start with an empty array + if len(lattice_neighbors) == 0: + return + + # Left edge + if i == 0: + pattern[:, 0] = lattice_neighbors[3] + if j == 0: + pattern[0, 0] = lattice_neighbors[0] + elif j == self.shape[1] - 1: + pattern[2, 0] = lattice_neighbors[5] + + # Bottom edge + if j == 0: + pattern[0, 1] = lattice_neighbors[1] + if i != 0: + pattern[0, 0] = lattice_neighbors[1] + if i != self.shape[0] - 1: + pattern[0, 2] = lattice_neighbors[1] + + # Right edge + if i == self.shape[0] - 1: + pattern[:, 2] = lattice_neighbors[4] + if j == 0: + pattern[0, 2] = lattice_neighbors[2] + elif j == self.shape[1] - 1: + pattern[2, 2] = lattice_neighbors[7] + + # Top edge + if j == self.shape[1] - 1: + pattern[2, 1] = lattice_neighbors[6] + if i != 0: + pattern[2, 0] = lattice_neighbors[6] + if i != self.shape[0] - 1: + pattern[2, 2] = lattice_neighbors[6] + + # Define an auxiliary function that returns a universe's neighbors + # among the universes inside the lattice + def find_lattice_neighbors(pattern, i, j): + + # Away from left edge + if i != 0: + if j > 0: + pattern[0, 0] = key(self.universes[j-1][i-1]) + pattern[1, 0] = key(self.universes[j][i-1]) + if j < self.shape[1] - 1: + pattern[2, 0] = key(self.universes[j+1][i-1]) + + # Away from bottom edge + if j != 0: + if i > 0: + pattern[0, 0] = key(self.universes[j-1][i-1]) + pattern[0, 1] = key(self.universes[j-1][i]) + if i < self.shape[0] - 1: + pattern[0, 2] = key(self.universes[j-1][i+1]) + + # Away from right edge + if i != self.shape[0] - 1: + if j > 0: + pattern[0, 2] = key(self.universes[j-1][i+1]) + pattern[1, 2] = key(self.universes[j][i+1]) + if j < self.shape[1] - 1: + pattern[2, 2] = key(self.universes[j+1][i+1]) + + # Away from top edge + if j != self.shape[1] - 1: + if i > 0: + pattern[2, 0] = key(self.universes[j+1][i-1]) + pattern[2, 1] = key(self.universes[j+1][i]) + if i < self.shape[0] - 1: + pattern[2, 2] = key(self.universes[j+1][i+1]) + + # Analyze lattice, find unique patterns in groups of universes + for j in range(self.shape[1]): + for i in range(self.shape[0]): + + # Skip universes to ignore + if self.universes[j][i] in universes_to_ignore: + continue + + # Create a neighborhood pattern based on the universe's + # neighbors in the grid, and lattice's neighbors at the edges + + # Degenerate discretization has all universes be different + if strategy == "degenerate": + patterns[(i, j)] = {'locations': [(i, j)]} + continue + + # Find neighbors among lattice's neighbors at the edges + find_edge_neighbors(pattern, i, j) + + # Find neighbors among the lattice's universes + find_lattice_neighbors(pattern, i, j) + + pattern[1, 1] = key(self.universes[j][i]) + + # Look for pattern in dictionary of patterns found + found = False + for known_pattern, pattern_data in patterns.items(): + + # Look at all rotations of pattern + for rot in range(4): + if not found and tuple(map(tuple, pattern)) ==\ + known_pattern: + found = True + + # Save location of the pattern in the lattice + pattern_data['locations'].append((i, j)) + + # Rotate pattern + pattern = np.rot90(pattern) + + # Look at transpose of pattern and its rotations + pattern = np.transpose(pattern) + for rot in range(4): + if not found and tuple(map(tuple, pattern)) ==\ + known_pattern: + found = True + + # Save location of the pattern in the lattice + pattern_data['locations'].append((i, j)) + + # Rotate pattern + pattern = np.rot90(pattern) + + # Transpose pattern back for the next search + pattern = np.transpose(pattern) + + # Create new pattern and add to the patterns dictionary + if not found: + patterns[tuple(map(tuple, pattern))] =\ + {'locations': [(i, j)]} + + # Discretize lattice + for pattern, pattern_data in patterns.items(): + + first_pos = pattern_data['locations'][0] + + # Create a clone of the universe, without cloning materials + new_universe = self.universes[first_pos[1]][first_pos[0]].clone( + clone_materials=False, clone_regions=False) + + # Replace only the materials in materials_to_clone + for material in materials_to_clone: + material_cloned = False + + for cell in new_universe.get_all_cells().values(): + + if cell.fill_type == 'material': + if cell.fill.id == material.id: + + # Only a single clone of each material is necessary + if not material_cloned: + material_clone = material.clone() + material_cloned = True + + cell.fill = material_clone + elif cell.fill_type == 'distribmat': + raise(ValueError, "Lattice discretization should not " + "be used with distributed materials") + elif len(cell.temperature) > 1 or len(cell.fill) > 1: + raise(ValueError, "Lattice discretization should not " + "be used with distributed cells") + + # Rebuild lattice from list of locations with this pattern + for index, location in enumerate(pattern_data['locations']): + self.universes[location[1]][location[0]] = new_universe + def create_xml_subelement(self, xml_element, memo=None): """Add the lattice xml representation to an incoming xml element diff --git a/openmc/universe.py b/openmc/universe.py index b64a877854..1d65397f86 100644 --- a/openmc/universe.py +++ b/openmc/universe.py @@ -482,12 +482,18 @@ class Universe(IDManagerMixin): return universes - def clone(self, memo=None): + def clone(self, clone_materials=True, clone_regions=True, memo=None): """Create a copy of this universe with a new unique ID, and clones all cells within this universe. Parameters ---------- + clone_materials : bool + Whether to create separates copies of the materials filling cells + contained in this universe. + clone_regions : bool + Whether to create separates copies of the regions bounding cells + contained in this universe. memo : dict or None A nested dictionary of previously cloned objects. This parameter is used internally and should not be specified by the user. @@ -510,7 +516,8 @@ class Universe(IDManagerMixin): # Clone all cells for the universe clone clone._cells = OrderedDict() for cell in self._cells.values(): - clone.add_cell(cell.clone(memo)) + clone.add_cell(cell.clone(clone_materials, clone_regions, + memo)) # Memoize the clone memo[self] = clone diff --git a/tests/unit_tests/test_cell.py b/tests/unit_tests/test_cell.py index d01d94a324..6ac4ae9c07 100644 --- a/tests/unit_tests/test_cell.py +++ b/tests/unit_tests/test_cell.py @@ -54,6 +54,18 @@ def test_clone(): assert c2.region != c.region assert c2.temperature == c.temperature + c3 = c.clone(clone_materials=False) + assert c3.id != c.id + assert c3.fill == c.fill + assert c3.region != c.region + assert c3.temperature == c.temperature + + c4 = c.clone(clone_regions=False) + assert c4.id != c.id + assert c4.fill != c.fill + assert c4.region == c.region + assert c4.temperature == c.temperature + def test_temperature(cell_with_lattice): # Make sure temperature propagates through universes diff --git a/tests/unit_tests/test_lattice.py b/tests/unit_tests/test_lattice.py index c4a3517769..31433af7dc 100644 --- a/tests/unit_tests/test_lattice.py +++ b/tests/unit_tests/test_lattice.py @@ -278,6 +278,16 @@ def test_clone(rlat2, hlat2, hlat3): assert hlat_clone.center == hlat3.center assert hlat_clone.pitch == hlat3.pitch + rlat_clone = rlat2.clone(clone_materials=False) + assert rlat_clone.get_all_materials() == rlat2.get_all_materials() + + rlat_clone = rlat2.clone(clone_materials=False, clone_regions=False) + for c1 in rlat_clone.cells: + for c2 in rlat2.cells: + if c1.fill == c2.fill: + print(c1.fill) + assert c1.region == c2.region + def test_repr(rlat2, rlat3, hlat2, hlat3): repr(rlat2) diff --git a/tests/unit_tests/test_lattice_discretization.py b/tests/unit_tests/test_lattice_discretization.py new file mode 100644 index 0000000000..b2766572d7 --- /dev/null +++ b/tests/unit_tests/test_lattice_discretization.py @@ -0,0 +1,91 @@ +from math import sqrt +import xml.etree.ElementTree as ET + +import openmc +import pytest + +from tests.unit_tests.test_lattice import zr, pincell1, pincell2, rlat2 + + +def test_discretization_degenerate(rlat2): + + rlat_clone = rlat2.clone() + rlat_clone.discretize() + + assert rlat_clone.get_universe((0, 0)).id !=\ + rlat_clone.get_universe((1, 0)).id + assert rlat_clone.get_universe((1, 0)).id !=\ + rlat_clone.get_universe((0, 1)).id + assert rlat_clone.get_universe((0, 1)).id !=\ + rlat_clone.get_universe((0, 0)).id + assert rlat_clone.get_universe((0, 2)).id !=\ + rlat_clone.get_universe((2, 0)).id + assert rlat_clone.get_universe((2, 1)).id !=\ + rlat_clone.get_universe((1, 2)).id + +def test_discretization_skip_universe(rlat2): + + rlat_clone = rlat2.clone() + rlat_clone.discretize( + universes_to_ignore=[rlat_clone.get_universe((0, 0))]) + + assert rlat_clone.get_universe((0, 0)) == rlat_clone.get_universe((0, 1)) + assert rlat_clone.get_universe((0, 1)) == rlat_clone.get_universe((2, 1)) + assert rlat_clone.get_universe((0, 1)) == rlat_clone.get_universe((1, 2)) + +def test_discretization_clone_only_some_materials(rlat2): + + rlat_clone = rlat2.clone() + fuel1 = next(iter(rlat_clone.get_universe((0, 0)).cells.values())).fill + rlat_clone.discretize(materials_to_clone=[fuel1]) + + assert next(reversed(rlat_clone.get_universe((0, 0)).cells.values())).fill\ + == next(reversed(rlat_clone.get_universe((1, 0)).cells.values())).fill + assert next(iter(rlat_clone.get_universe((0, 0)).cells.values())).fill\ + != next(iter(rlat_clone.get_universe((1, 0)).cells.values())).fill + +def test_discretization_lns(rlat2): + + rlat_clone = rlat2.clone() + rlat_clone.discretize(strategy="lns") + + assert rlat_clone.get_universe((0, 2)) == rlat_clone.get_universe((2, 2)) + assert rlat_clone.get_universe((0, 2)) != rlat_clone.get_universe((0, 0)) + assert rlat_clone.get_universe((0, 2)) != rlat_clone.get_universe((0, 1)) + assert rlat_clone.get_universe((0, 2)) != rlat_clone.get_universe((1, 0)) + assert rlat_clone.get_universe((0, 2)) != rlat_clone.get_universe((1, 1)) + assert rlat_clone.get_universe((0, 2)) != rlat_clone.get_universe((1, 2)) + assert rlat_clone.get_universe((0, 2)) != rlat_clone.get_universe((2, 0)) + assert rlat_clone.get_universe((0, 2)) != rlat_clone.get_universe((2, 1)) + +def test_discretization_lns_using_names(rlat2): + + rlat_clone = rlat2.clone() + rlat_clone.discretize(strategy="lns", key=lambda univ: univ.name) + + assert rlat_clone.get_universe((0, 2)) == rlat_clone.get_universe((2, 2)) + + rlat_clone = rlat2.clone() + rlat_clone.get_universe((0, 1)).name="u1" + rlat_clone.get_universe((1, 0)).name="u1" + rlat_clone.discretize(strategy="lns", key=lambda univ: univ.name) + + assert rlat_clone.get_universe((0, 2)) == rlat_clone.get_universe((2, 0)) + assert rlat_clone.get_universe((2, 2)) == rlat_clone.get_universe((0, 0)) + assert rlat_clone.get_universe((0, 2)) == rlat_clone.get_universe((2, 2)) + assert rlat_clone.get_universe((1, 2)) == rlat_clone.get_universe((1, 0)) + assert rlat_clone.get_universe((2, 1)) == rlat_clone.get_universe((0, 1)) + assert rlat_clone.get_universe((1, 2)) == rlat_clone.get_universe((0, 1)) + assert rlat_clone.get_universe((0, 0)) != rlat_clone.get_universe((0, 1)) + +def test_discretization_lns_with_neighbor_list(rlat2): + + rlat_clone = rlat2.clone() + u1 = rlat_clone.get_universe((1, 0)) + u2 = rlat_clone.get_universe((0, 0)) + rlat_clone.discretize(strategy="lns", + lattice_neighbors=[u1,u1,u2,u1,u2,u2,u1,u2]) + + assert rlat_clone.get_universe((0, 1)) == rlat_clone.get_universe((0, 0)) + assert rlat_clone.get_universe((1, 1)) != rlat_clone.get_universe((0, 0)) + diff --git a/tests/unit_tests/test_universe.py b/tests/unit_tests/test_universe.py index 674fc5c185..ab3fd6e33a 100644 --- a/tests/unit_tests/test_universe.py +++ b/tests/unit_tests/test_universe.py @@ -99,6 +99,28 @@ def test_get_all_universes(): assert not (univs ^ {u1, u2}) +def test_clone(): + + c1 = openmc.Cell(cell_id=1) + c1.region = -openmc.ZCylinder(r=1.0) + c2 = openmc.Cell(cell_id=2) + c2.fill = openmc.Material() + c3 = openmc.Cell() + u1 = openmc.Universe(name='cool', cells=(c1, c2, c3)) + + u2 = u1.clone() + assert u2.name == u1.name + assert u2.cells != u1.cells + assert u2.get_all_materials() != u1.get_all_materials() + + u2 = u1.clone(clone_materials=False) + assert u2.get_all_materials() == u1.get_all_materials() + + u3 = u1.clone(clone_regions=False) + assert next(iter(u3.cells.values())).region ==\ + next(iter(u1.cells.values())).region + + def test_create_xml(cell_with_lattice): cells = [openmc.Cell() for i in range(5)] u = openmc.Universe(cells=cells)