diff --git a/docs/source/_images/stack_lat.svg b/docs/source/_images/stack_lat.svg new file mode 100644 index 00000000000..df9c66c6c81 --- /dev/null +++ b/docs/source/_images/stack_lat.svg @@ -0,0 +1,253 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/docs/source/methods/geometry.rst b/docs/source/methods/geometry.rst index b282ffdee37..a1cad38436e 100644 --- a/docs/source/methods/geometry.rst +++ b/docs/source/methods/geometry.rst @@ -549,6 +549,45 @@ volumetric units of the lattice will be referred to as "tiles". Tiles are identified by thier indices, and the process of discovering which tile contains the particle is referred to as "indexing". +Stack Lattice Indexing +---------------------- + +Indices are assigned to tiles in a stack lattice based on the tile's +position along the axis of orientation (the :math:`x`, :math:`y`, or +:math:`z` axis). Figure :num:`fig-stack-lat` maps the indices for a +stack lattice. The index, (1), maps to the bottom tile. (5) maps to +the top tile. + + +.. _fig-stack-lat: + +.. figure:: ../_images/stack_lat.svg + :align: center + :figclass: align-center + :width: 400px + + Stack lattice tile indices. + +In general, a lattice tile is specified by the single index, +:math:`(i)`. If a particle's current coordinates are :math:`(x, y, z)` +and the lattice is oriented in the :math:`z` direction, then the index +of a tile in a uniform stack lattice can be determined from this formula: + +.. math:: + :label: stack_indexing + + i = \left \lfloor \frac{z - z_0}{p} \right \rfloor + +where :math:`(z_0)` is the coordinate to the base of the lattice, and +:math:`p` is the pitch along the :math:`z` axis. Similar equations govern +uniform stack lattices aligned with the :math:`x` and :math:`y` axes. + + +The index of a tile in a nonuniform stack lattice (i.e. a stack lattice where +tiles have differing dimensions) does not have a closed-form equation. It can +only be found numerically. + + Rectilinear Lattice Indexing ---------------------------- diff --git a/docs/source/pythonapi/base.rst b/docs/source/pythonapi/base.rst index b471d12a86f..ffd4967f03d 100644 --- a/docs/source/pythonapi/base.rst +++ b/docs/source/pythonapi/base.rst @@ -92,6 +92,7 @@ Building geometry openmc.DAGMCUniverse openmc.RectLattice openmc.HexLattice + openmc.StackLattice openmc.Geometry Many of the above classes are derived from several abstract classes: diff --git a/docs/source/usersguide/geometry.rst b/docs/source/usersguide/geometry.rst index 2382609a48c..b08cfc5c60f 100644 --- a/docs/source/usersguide/geometry.rst +++ b/docs/source/usersguide/geometry.rst @@ -287,6 +287,9 @@ regular pattern such as a rectangular or hexagonal lattice. In such a case, it would be cumbersome to have to define the boundaries of each of the cells to be filled with a universe. OpenMC provides a means to define lattice structures through the :class:`openmc.RectLattice` and :class:`openmc.HexLattice` classes. +Sometimes a model will have a repeated structure in only one dimension. OpenMC +provides a means to define such structures through the +:class:`openmc.StackLattice` class. Rectangular Lattices -------------------- @@ -412,6 +415,28 @@ code would work:: If you need to create a hexagonal boundary (composed of six planar surfaces) for a hexagonal lattice, :func:`openmc.model.hexagonal_prism` can be used. +Stack Lattices +-------------- +A stack lattice defines a one-dimensional array of universes that are filled into stacked layers (lattice elements). These layers can be of uniform or nonuniform thickness. Creating a stack lattice is similar to creating a rectangular lattice with a few differences: + +The position of the central axis and coordinate of the base level of the stack +in must be specified (:attr:`StackLattice.central_axis` and +:attr:`StackLattice.base_coordinate`, respectively). Additionally, +:attr:`universes` must be set before :attr:`pitch`. For nonuniform lattices, the :attr:`is_uniform` attribute must be set to ``False`` before setting :attr:`pitch`. + +For a uniform stack lattice, a single value for the pitch should be specified. +It should not be a list. For a nonuniform stack lattice, an iterable of values +should be given. The length of this iterable should match the length of +:attr:`StackLattice.universes`. + +The lattice orientation is by default in the :math:`z`-direction, but this can +be changed using the :attr:`StackLattice.orientation` attribute to :attr:`x` or +:attr:`y`. :attr:`StackLattice.central_axis` assumes that the coordinates are +in the :math:`(x,y)`, :math:`(x,z)`, or :math:`(y,z)` basis for +:attr:`orientation` of :attr:`z`, :attr:`y`, or :attr:`x`, respectively. + + + .. _usersguide_geom_export: -------------------------- diff --git a/include/openmc/hdf5_interface.h b/include/openmc/hdf5_interface.h index e9d9b4f9682..1916edce03f 100644 --- a/include/openmc/hdf5_interface.h +++ b/include/openmc/hdf5_interface.h @@ -52,6 +52,7 @@ inline hid_t create_group(hid_t parent_id, const std::stringstream& name) hid_t file_open(const std::string& filename, char mode, bool parallel = false); hid_t open_group(hid_t group_id, const std::string& name); +void write_bool(hid_t group_id, const char* name, bool b, bool indep); void write_string( hid_t group_id, const char* name, const std::string& buffer, bool indep); @@ -109,6 +110,8 @@ void write_double(hid_t group_id, int ndim, const hsize_t* dims, const char* name, const double* buffer, bool indep); void write_int(hid_t group_id, int ndim, const hsize_t* dims, const char* name, const int* buffer, bool indep); +void write_bool(hid_t group_id, int ndim, const hsize_t* dims, const char* name, + const hbool_t* buffer, bool indep); void write_llong(hid_t group_id, int ndim, const hsize_t* dims, const char* name, const long long* buffer, bool indep); void write_string(hid_t group_id, int ndim, const hsize_t* dims, size_t slen, diff --git a/include/openmc/lattice.h b/include/openmc/lattice.h index a2411467df1..27ff727ecec 100644 --- a/include/openmc/lattice.h +++ b/include/openmc/lattice.h @@ -22,7 +22,7 @@ namespace openmc { constexpr int32_t NO_OUTER_UNIVERSE {-1}; -enum class LatticeType { rect, hex }; +enum class LatticeType { rect, hex, stack }; //============================================================================== // Global variables @@ -285,6 +285,51 @@ class HexLattice : public Lattice { array pitch_; //!< Lattice tile width and height }; +//============================================================================== + +class StackLattice : public Lattice { +public: + explicit StackLattice(pugi::xml_node lat_node); + + int32_t const& operator[](array const& i_xyz); + + bool are_valid_indices(array const& i_xyz) const; + + std::pair> distance( + Position r, Direction u, const array& i_xyz) const; + + void get_indices(Position r, Direction u, array& result) const; + + int get_flat_index(const array& i_xyz) const; + + Position get_local_position(Position r, const array& i_xyz) const; + + int32_t& offset(int map, array const& i_xyz); + + int32_t offset(int map, int indx) const; + + std::string index_to_string(int indx) const; + + void to_hdf5_inner(hid_t group_id) const; + +private: + enum class Orientation { + z, //!< Central axis of lattice parallel to z-axis + y, //!< Central axis of lattice parallel to y-axis + x //!< Central axis of lattice parallel to x-axis + }; + + bool is_uniform_; //!< Mark if lattice is uniform or not + + int orientation_idx_; //!< Index to select x,y, or z based on orientation + int n_layers_; //!< Number of layer positions + Orientation orientation_; //!< Orientation of lattice + Position central_axis_; //!< Axial center of lattice + double base_coordinate_; //!< Coordinate of base layer of lattice + vector pitch_; //!< Lattice layer width. + vector layer_boundaries_; //!< Coordinates of lattice tile boundaries +}; + //============================================================================== // Non-member functions //============================================================================== diff --git a/openmc/geometry.py b/openmc/geometry.py index 58a33c1b121..cb5f3d934aa 100644 --- a/openmc/geometry.py +++ b/openmc/geometry.py @@ -110,7 +110,8 @@ def export_to_xml(self, path='geometry.xml', remove_surfs=False): p /= 'geometry.xml' # Write the XML Tree to the geometry.xml file - xml.reorder_attributes(root_element) # TODO: Remove when support is Python 3.8+ + # TODO: Remove when support is Python 3.8+ + xml.reorder_attributes(root_element) tree = ET.ElementTree(root_element) tree.write(str(p), xml_declaration=True, encoding='utf-8') @@ -134,6 +135,7 @@ def from_xml(cls, path='geometry.xml', materials=None): """ # Helper function for keeping a cache of Universe instances universes = {} + def get_universe(univ_id): if univ_id not in universes: univ = openmc.Universe(univ_id) @@ -191,6 +193,14 @@ def get_universe(univ_id): for u in ring: child_of[u].append(lat) + for elem in root.findall('stack_lattice'): + lat = openmc.StackLattice.from_xml_element(elem, get_universe) + universes[lat.id] = lat + if lat.outer is not None: + child_of[lat.outer].append(lat) + for u in lat.universes.ravel(): + child_of[u].append(lat) + # Create dictionary to easily look up materials if materials is None: filename = Path(path).parent / 'materials.xml' @@ -199,7 +209,8 @@ def get_universe(univ_id): mats['void'] = None for elem in root.findall('cell'): - c = openmc.Cell.from_xml_element(elem, surfaces, mats, get_universe) + c = openmc.Cell.from_xml_element( + elem, surfaces, mats, get_universe) if c.fill_type in ('universe', 'lattice'): child_of[c.fill].append(c) @@ -258,7 +269,7 @@ def get_instances(self, paths): for p in path_list: # Extract the cell id from the path last_index = p.rfind('>') - last_path = p[last_index+1:] + last_path = p[last_index + 1:] uid = int(last_path[1:]) # Get corresponding cell/material @@ -418,7 +429,12 @@ def get_redundant_surfaces(self): for keep, *redundant in tally.values() for replace in redundant} - def _get_domains_by_name(self, name, case_sensitive, matching, domain_type): + def _get_domains_by_name( + self, + name, + case_sensitive, + matching, + domain_type): if not case_sensitive: name = name.lower() @@ -435,7 +451,11 @@ def _get_domains_by_name(self, name, case_sensitive, matching, domain_type): domains.sort(key=lambda x: x.id) return domains - def get_materials_by_name(self, name, case_sensitive=False, matching=False): + def get_materials_by_name( + self, + name, + case_sensitive=False, + matching=False): """Return a list of materials with matching names. Parameters @@ -454,7 +474,8 @@ def get_materials_by_name(self, name, case_sensitive=False, matching=False): Materials matching the queried name """ - return self._get_domains_by_name(name, case_sensitive, matching, 'material') + return self._get_domains_by_name( + name, case_sensitive, matching, 'material') def get_cells_by_name(self, name, case_sensitive=False, matching=False): """Return a list of cells with matching names. @@ -475,9 +496,14 @@ def get_cells_by_name(self, name, case_sensitive=False, matching=False): Cells matching the queried name """ - return self._get_domains_by_name(name, case_sensitive, matching, 'cell') - - def get_cells_by_fill_name(self, name, case_sensitive=False, matching=False): + return self._get_domains_by_name( + name, case_sensitive, matching, 'cell') + + def get_cells_by_fill_name( + self, + name, + case_sensitive=False, + matching=False): """Return a list of cells with fills with matching names. Parameters @@ -522,7 +548,11 @@ def get_cells_by_fill_name(self, name, case_sensitive=False, matching=False): return sorted(cells, key=lambda x: x.id) - def get_universes_by_name(self, name, case_sensitive=False, matching=False): + def get_universes_by_name( + self, + name, + case_sensitive=False, + matching=False): """Return a list of universes with matching names. Parameters @@ -541,7 +571,8 @@ def get_universes_by_name(self, name, case_sensitive=False, matching=False): Universes matching the queried name """ - return self._get_domains_by_name(name, case_sensitive, matching, 'universe') + return self._get_domains_by_name( + name, case_sensitive, matching, 'universe') def get_lattices_by_name(self, name, case_sensitive=False, matching=False): """Return a list of lattices with matching names. @@ -562,7 +593,8 @@ def get_lattices_by_name(self, name, case_sensitive=False, matching=False): Lattices matching the queried name """ - return self._get_domains_by_name(name, case_sensitive, matching, 'lattice') + return self._get_domains_by_name( + name, case_sensitive, matching, 'lattice') def remove_redundant_surfaces(self): """Remove redundant surfaces from the geometry""" diff --git a/openmc/lattice.py b/openmc/lattice.py index 464daaf7e51..60c87159f3a 100644 --- a/openmc/lattice.py +++ b/openmc/lattice.py @@ -105,6 +105,8 @@ def from_hdf5(group, universes): return openmc.RectLattice.from_hdf5(group, universes) elif lattice_type == 'hexagonal': return openmc.HexLattice.from_hdf5(group, universes) + elif lattice_type == 'stack': + return openmc.StackLattice.from_hdf5(group, universes) else: raise ValueError(f'Unknown lattice type: {lattice_type}') @@ -121,15 +123,19 @@ def get_unique_universes(self): univs = OrderedDict() for k in range(len(self._universes)): - for j in range(len(self._universes[k])): - if isinstance(self._universes[k][j], openmc.Universe): - u = self._universes[k][j] - univs[u._id] = u - else: - for i in range(len(self._universes[k][j])): - u = self._universes[k][j][i] - assert isinstance(u, openmc.Universe) + if isinstance(self._universes[k], openmc.Universe): + u = self._universes[k] + univs[u._id] = u + else: + for j in range(len(self._universes[k])): + if isinstance(self._universes[k][j], openmc.Universe): + u = self._universes[k][j] univs[u._id] = u + else: + for i in range(len(self._universes[k][j])): + u = self._universes[k][j][i] + assert isinstance(u, openmc.Universe) + univs[u._id] = u if self.outer is not None: univs[self.outer._id] = self.outer @@ -242,7 +248,8 @@ def get_universe(self, idx): hexagonal lattices, they are given in the :math:`x,\alpha` or :math:`x,\alpha,z` coordinate systems for "y" orientations and :math:`\alpha,y` or :math:`\alpha,y,z` coordinate systems for "x" - orientations. + orientations. For stack lattices, they are given in the coordinate + corresponding to the orientation axis. Returns ------- @@ -251,6 +258,8 @@ def get_universe(self, idx): """ idx_u = self.get_universe_index(idx) + if self.ndim == 1: + return self.universes[idx_u] if self.ndim == 2: return self.universes[idx_u[0]][idx_u[1]] else: @@ -314,18 +323,22 @@ def clone(self, clone_materials=True, clone_regions=True, memo=None): if self.outer is not None: clone.outer = self.outer.clone(clone_materials, clone_regions, - memo) + memo) # Assign universe clones to the lattice clone for i in self.indices: if isinstance(self, RectLattice): clone.universes[i] = self.universes[i].clone( - clone_materials, clone_regions, memo) + clone_materials, clone_regions, memo) else: - if self.ndim == 2: + if self.ndim == 1: + clone.universes[i] = \ + self.universes[i].clone(clone_materials, + clone_regions, memo) + elif self.ndim == 2: clone.universes[i[0]][i[1]] = \ self.universes[i[0]][i[1]].clone(clone_materials, - clone_regions, memo) + clone_regions, memo) else: clone.universes[i[0]][i[1]][i[2]] = \ self.universes[i[0]][i[1]][i[2]].clone( @@ -488,8 +501,12 @@ def pitch(self, pitch): @Lattice.universes.setter def universes(self, universes): - cv.check_iterable_type('lattice universes', universes, openmc.UniverseBase, - min_depth=2, max_depth=3) + cv.check_iterable_type( + 'lattice universes', + universes, + openmc.UniverseBase, + min_depth=2, + max_depth=3) self._universes = np.asarray(universes) def find_element(self, point): @@ -509,12 +526,12 @@ def find_element(self, point): element coordinate system """ - ix = floor((point[0] - self.lower_left[0])/self.pitch[0]) - iy = floor((point[1] - self.lower_left[1])/self.pitch[1]) + ix = floor((point[0] - self.lower_left[0]) / self.pitch[0]) + iy = floor((point[1] - self.lower_left[1]) / self.pitch[1]) if self.ndim == 2: idx = (ix, iy) else: - iz = floor((point[2] - self.lower_left[2])/self.pitch[2]) + iz = floor((point[2] - self.lower_left[2]) / self.pitch[2]) idx = (ix, iy, iz) return idx, self.get_local_coordinates(point, idx) @@ -536,12 +553,13 @@ def get_local_coordinates(self, point, idx): system """ - x = point[0] - (self.lower_left[0] + (idx[0] + 0.5)*self.pitch[0]) - y = point[1] - (self.lower_left[1] + (idx[1] + 0.5)*self.pitch[1]) + x = point[0] - (self.lower_left[0] + (idx[0] + 0.5) * self.pitch[0]) + y = point[1] - (self.lower_left[1] + (idx[1] + 0.5) * self.pitch[1]) if self.ndim == 2: z = point[2] else: - z = point[2] - (self.lower_left[2] + (idx[2] + 0.5)*self.pitch[2]) + z = point[2] - (self.lower_left[2] + + (idx[2] + 0.5) * self.pitch[2]) return (x, y, z) def get_universe_index(self, idx): @@ -703,34 +721,34 @@ 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]) + 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]) + 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]) + 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]) + 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]) + 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]) + 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]) + 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]) + 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]): @@ -763,7 +781,7 @@ def find_lattice_neighbors(pattern, i, j): # Look at all rotations of pattern for rot in range(4): if not found and tuple(map(tuple, pattern)) ==\ - known_pattern: + known_pattern: found = True # Save location of the pattern in the lattice @@ -776,7 +794,7 @@ def find_lattice_neighbors(pattern, i, j): pattern = np.transpose(pattern) for rot in range(4): if not found and tuple(map(tuple, pattern)) ==\ - known_pattern: + known_pattern: found = True # Save location of the pattern in the lattice @@ -791,7 +809,7 @@ def find_lattice_neighbors(pattern, i, j): # Create new pattern and add to the patterns dictionary if not found: patterns[tuple(map(tuple, pattern))] =\ - {'locations': [(i, j)]} + {'locations': [(i, j)]} # Discretize lattice for pattern, pattern_data in patterns.items(): @@ -800,7 +818,7 @@ def find_lattice_neighbors(pattern, i, j): # 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) + clone_materials=False, clone_regions=False) # Replace only the materials in materials_to_clone for material in materials_to_clone: @@ -854,7 +872,8 @@ def create_xml_subelement(self, xml_element, memo=None): # Make sure universes have been assigned if self.universes is None: - raise ValueError(f"Lattice {self.id} does not have universes assigned.") + raise ValueError( + f"Lattice {self.id} does not have universes assigned.") lattice_subelement = ET.Element("lattice") lattice_subelement.set("id", str(self._id)) @@ -1096,8 +1115,10 @@ def __repr__(self): string += '{0: <16}{1}{2}\n'.format('\tName', '=\t', self._name) string += '{0: <16}{1}{2}\n'.format('\tOrientation', '=\t', self._orientation) - string += '{0: <16}{1}{2}\n'.format('\t# Rings', '=\t', self._num_rings) - string += '{0: <16}{1}{2}\n'.format('\t# Axial', '=\t', self._num_axial) + string += '{0: <16}{1}{2}\n'.format('\t# Rings', + '=\t', self._num_rings) + string += '{0: <16}{1}{2}\n'.format('\t# Axial', + '=\t', self._num_axial) string += '{0: <16}{1}{2}\n'.format('\tCenter', '=\t', self._center) string += '{0: <16}{1}{2}\n'.format('\tPitch', '=\t', self._pitch) @@ -1140,11 +1161,11 @@ def center(self): def indices(self): if self.num_axial is None: return [(r, i) for r in range(self.num_rings) - for i in range(max(6*(self.num_rings - 1 - r), 1))] + for i in range(max(6 * (self.num_rings - 1 - r), 1))] else: return [(z, r, i) for z in range(self.num_axial) for r in range(self.num_rings) - for i in range(max(6*(self.num_rings - 1 - r), 1))] + for i in range(max(6 * (self.num_rings - 1 - r), 1))] @property def _natural_indices(self): @@ -1236,13 +1257,13 @@ def universes(self, universes): raise ValueError(msg) # Check the outer rings. - for r in range(self._num_rings-1): - if len(axial_slice[r]) != 6*(self._num_rings - 1 - r): + for r in range(self._num_rings - 1): + if len(axial_slice[r]) != 6 * (self._num_rings - 1 - r): msg = 'HexLattice ID={0:d} has the wrong number of ' \ 'elements in ring number {1:d} (counting from the '\ 'outermost ring). This ring should have {2:d} ' \ 'elements.'.format(self._id, r, - 6*(self._num_rings - 1 - r)) + 6 * (self._num_rings - 1 - r)) raise ValueError(msg) else: @@ -1255,13 +1276,13 @@ def universes(self, universes): raise ValueError(msg) # Check the outer rings. - for r in range(self._num_rings-1): - if len(axial_slice[r]) != 6*(self._num_rings - 1 - r): + for r in range(self._num_rings - 1): + if len(axial_slice[r]) != 6 * (self._num_rings - 1 - r): msg = 'HexLattice ID={0:d} has the wrong number of ' \ 'elements in ring number {1:d} (counting from the '\ 'outermost ring). This ring should have {2:d} ' \ 'elements.'.format(self._id, r, - 6*(self._num_rings - 1 - r)) + 6 * (self._num_rings - 1 - r)) raise ValueError(msg) def find_element(self, point): @@ -1289,15 +1310,15 @@ def find_element(self, point): iz = 1 else: z = point[2] - self.center[2] - iz = floor(z/self.pitch[1] + 0.5*self.num_axial) + iz = floor(z / self.pitch[1] + 0.5 * self.num_axial) if self._orientation == 'x': - alpha = y - x*sqrt(3.) - i1 = floor(-alpha/(sqrt(3.0) * self.pitch[0])) - i2 = floor(y/(sqrt(0.75) * self.pitch[0])) + alpha = y - x * sqrt(3.) + i1 = floor(-alpha / (sqrt(3.0) * self.pitch[0])) + i2 = floor(y / (sqrt(0.75) * self.pitch[0])) else: - alpha = y - x/sqrt(3.) - i1 = floor(x/(sqrt(0.75) * self.pitch[0])) - i2 = floor(alpha/self.pitch[0]) + alpha = y - x / sqrt(3.) + i1 = floor(x / (sqrt(0.75) * self.pitch[0])) + i2 = floor(alpha / self.pitch[0]) # Check four lattice elements to see which one is closest based on local # coordinates indices = [(i1, i2, iz), (i1 + 1, i2, iz), (i1, i2 + 1, iz), @@ -1333,17 +1354,21 @@ def get_local_coordinates(self, point, idx): """ if self._orientation == 'x': - x = point[0] - (self.center[0] + (idx[0] + 0.5*idx[1])*self.pitch[0]) - y = point[1] - (self.center[1] + sqrt(0.75)*self.pitch[0]*idx[1]) + x = point[0] - (self.center[0] + (idx[0] + + 0.5 * idx[1]) * self.pitch[0]) + y = point[1] - (self.center[1] + sqrt(0.75) + * self.pitch[0] * idx[1]) else: - x = point[0] - (self.center[0] + sqrt(0.75)*self.pitch[0]*idx[0]) - y = point[1] - (self.center[1] + (0.5*idx[0] + idx[1])*self.pitch[0]) + x = point[0] - (self.center[0] + sqrt(0.75) + * self.pitch[0] * idx[0]) + y = point[1] - (self.center[1] + + (0.5 * idx[0] + idx[1]) * self.pitch[0]) if self._num_axial is None: z = point[2] else: - z = point[2] - (self.center[2] + (idx[2] + 0.5 - 0.5*self.num_axial) * - self.pitch[1]) + z = point[2] - (self.center[2] + (idx[2] + 0.5 - + 0.5 * self.num_axial) * self.pitch[1]) return (x, y, z) def get_universe_index(self, idx): @@ -1371,21 +1396,22 @@ def get_universe_index(self, idx): z = -a - x g = max(abs(x), abs(a), abs(z)) - # Next we use a clever method to figure out where along the ring we are. + # Next we use a clever method to figure out where along the ring we + # are. i_ring = self._num_rings - 1 - g if x >= 0: if a >= 0: i_within = x else: - i_within = 2*g + z + i_within = 2 * g + z else: if a <= 0: - i_within = 3*g - x + i_within = 3 * g - x else: - i_within = 5*g - z + i_within = 5 * g - z if self._orientation == 'x' and g > 0: - i_within = (i_within + 5*g) % (6*g) + i_within = (i_within + 5 * g) % (6 * g) if self.num_axial is None: return (i_ring, i_within) @@ -1453,7 +1479,8 @@ def create_xml_subelement(self, xml_element, memo=None): # Export the Lattice nested Universe IDs. if self.universes is None: - raise ValueError(f"Lattice {self.id} does not have universes assigned.") + raise ValueError( + f"Lattice {self.id} does not have universes assigned.") # 3D Lattices if self._num_axial is not None: @@ -1464,8 +1491,8 @@ def create_xml_subelement(self, xml_element, memo=None): universe.create_xml_subelement(xml_element, memo) # Initialize the remaining universes. - for r in range(self._num_rings-1): - for theta in range(6*(self._num_rings - 1 - r)): + for r in range(self._num_rings - 1): + for theta in range(6 * (self._num_rings - 1 - r)): universe = self._universes[z][r][theta] universe.create_xml_subelement(xml_element, memo) @@ -1483,7 +1510,7 @@ def create_xml_subelement(self, xml_element, memo=None): # Initialize the remaining universes. for r in range(self._num_rings - 1): - for theta in range(6*(self._num_rings - 1 - r)): + for theta in range(6 * (self._num_rings - 1 - r)): universe = self._universes[r][theta] universe.create_xml_subelement(xml_element, memo) @@ -1529,7 +1556,7 @@ def from_xml_element(cls, elem, get_universe): lat._num_axial = n_axial = int(get_text(elem, 'n_axial', 1)) # Create empty nested lists for one axial level - univs = [[None for _ in range(max(6*(n_rings - 1 - r), 1))] + univs = [[None for _ in range(max(6 * (n_rings - 1 - r), 1))] for r in range(n_rings)] if n_axial > 1: univs = [deepcopy(univs) for i in range(n_axial)] @@ -1620,11 +1647,11 @@ def _repr_axial_slice_x(self, universes): largest_id = max([max([univ._id for univ in ring]) for ring in universes]) n_digits = len(str(largest_id)) - pad = ' '*n_digits + pad = ' ' * n_digits id_form = '{: ^' + str(n_digits) + 'd}' # Initialize the list for each row. - rows = [[] for i in range(2*self._num_rings - 1)] + rows = [[] for i in range(2 * self._num_rings - 1)] middle = self._num_rings - 1 # Start with the degenerate first ring. @@ -1701,8 +1728,8 @@ def _repr_axial_slice_x(self, universes): # Pad the beginning of the rows so they line up properly. for y in range(self._num_rings - 1): - rows[y] = (self._num_rings - 1 - y)*pad + rows[y] - rows[-1 - y] = (self._num_rings - 1 - y)*pad + rows[-1 - y] + rows[y] = (self._num_rings - 1 - y) * pad + rows[y] + rows[-1 - y] = (self._num_rings - 1 - y) * pad + rows[-1 - y] # Join the rows together and return the string. universe_ids = '\n'.join(rows) @@ -1722,11 +1749,11 @@ def _repr_axial_slice_y(self, universes): largest_id = max([max([univ._id for univ in ring]) for ring in universes]) n_digits = len(str(largest_id)) - pad = ' '*n_digits + pad = ' ' * n_digits id_form = '{: ^' + str(n_digits) + 'd}' # Initialize the list for each row. - rows = [[] for i in range(1 + 4 * (self._num_rings-1))] + rows = [[] for i in range(1 + 4 * (self._num_rings - 1))] middle = 2 * (self._num_rings - 1) # Start with the degenerate first ring. @@ -1738,7 +1765,7 @@ def _repr_axial_slice_y(self, universes): # r_prime increments down while r increments up. r_prime = self._num_rings - 1 - r theta = 0 - y = middle + 2*r + y = middle + 2 * r # Climb down the top-right. for i in range(r): @@ -1805,8 +1832,8 @@ def _repr_axial_slice_y(self, universes): # Pad the beginning of the rows so they line up properly. for y in range(self._num_rings - 1): - rows[y] = (self._num_rings - 1 - y)*pad + rows[y] - rows[-1 - y] = (self._num_rings - 1 - y)*pad + rows[-1 - y] + rows[y] = (self._num_rings - 1 - y) * pad + rows[y] + rows[-1 - y] = (self._num_rings - 1 - y) * pad + rows[-1 - y] for y in range(self._num_rings % 2, self._num_rings, 2): rows[middle + y] = pad + rows[middle + y] @@ -1849,14 +1876,14 @@ def _show_indices_y(num_rings): # Find the largest string and count the number of digits so we can # properly pad the output string later - largest_index = 6*(num_rings - 1) + largest_index = 6 * (num_rings - 1) n_digits_index = len(str(largest_index)) n_digits_ring = len(str(num_rings - 1)) str_form = '({{:{}}},{{:{}}})'.format(n_digits_ring, n_digits_index) - pad = ' '*(n_digits_index + n_digits_ring + 3) + pad = ' ' * (n_digits_index + n_digits_ring + 3) # Initialize the list for each row. - rows = [[] for i in range(1 + 4 * (num_rings-1))] + rows = [[] for i in range(1 + 4 * (num_rings - 1))] middle = 2 * (num_rings - 1) # Start with the degenerate first ring. @@ -1867,7 +1894,7 @@ def _show_indices_y(num_rings): # r_prime increments down while r increments up. r_prime = num_rings - 1 - r theta = 0 - y = middle + 2*r + y = middle + 2 * r for i in range(r): # Climb down the top-right. @@ -1910,8 +1937,8 @@ def _show_indices_y(num_rings): # Pad the beginning of the rows so they line up properly. for y in range(num_rings - 1): - rows[y] = (num_rings - 1 - y)*pad + rows[y] - rows[-1 - y] = (num_rings - 1 - y)*pad + rows[-1 - y] + rows[y] = (num_rings - 1 - y) * pad + rows[y] + rows[-1 - y] = (num_rings - 1 - y) * pad + rows[-1 - y] for y in range(num_rings % 2, num_rings, 2): rows[middle + y] = pad + rows[middle + y] @@ -1954,14 +1981,14 @@ def _show_indices_x(num_rings): # Find the largest string and count the number of digits so we can # properly pad the output string later - largest_index = 6*(num_rings - 1) + largest_index = 6 * (num_rings - 1) n_digits_index = len(str(largest_index)) n_digits_ring = len(str(num_rings - 1)) str_form = '({{:{}}},{{:{}}})'.format(n_digits_ring, n_digits_index) - pad = ' '*(n_digits_index + n_digits_ring + 3) + pad = ' ' * (n_digits_index + n_digits_ring + 3) # Initialize the list for each row. - rows = [[] for i in range(2*num_rings - 1)] + rows = [[] for i in range(2 * num_rings - 1)] middle = num_rings - 1 # Start with the degenerate first ring. @@ -2013,8 +2040,8 @@ def _show_indices_x(num_rings): # Pad the beginning of the rows so they line up properly. for y in range(num_rings - 1): - rows[y] = (num_rings - 1 - y)*pad + rows[y] - rows[-1 - y] = (num_rings - 1 - y)*pad + rows[-1 - y] + rows[y] = (num_rings - 1 - y) * pad + rows[y] + rows[-1 - y] = (num_rings - 1 - y) * pad + rows[-1 - y] # Join the rows together and return the string. return '\n\n'.join(rows) @@ -2091,7 +2118,7 @@ def from_hdf5(cls, group, universes): # Add a list for this axial level. uarray.append([]) x = n_rings - 1 - a = 2*n_rings - 2 + a = 2 * n_rings - 2 for r in range(n_rings - 1, 0, -1): # Add a list for this ring. uarray[-1].append([]) @@ -2133,7 +2160,7 @@ def from_hdf5(cls, group, universes): # Convert the ids into Universe objects. uarray[-1][-1] = [universes[u_id] - for u_id in uarray[-1][-1]] + for u_id in uarray[-1][-1]] # Handle the degenerate center ring separately. u_id = universe_ids[z, a, x] @@ -2147,7 +2174,7 @@ def from_hdf5(cls, group, universes): for z in range(n_axial): # Add a list for this axial level. uarray.append([]) - a = 2*n_rings - 2 + a = 2 * n_rings - 2 y = n_rings - 1 for r in range(n_rings - 1, 0, -1): # Add a list for this ring. @@ -2190,7 +2217,7 @@ def from_hdf5(cls, group, universes): # Convert the ids into Universe objects. uarray[-1][-1] = [universes[u_id] - for u_id in uarray[-1][-1]] + for u_id in uarray[-1][-1]] # Handle the degenerate center ring separately. u_id = universe_ids[z, y, a] @@ -2205,3 +2232,497 @@ def from_hdf5(cls, group, universes): lattice.universes = uarray[0] return lattice + + +class StackLattice(Lattice): + """A lattice consisting of universes stacked in layers in one dimension. + + To completely define a stack lattice, the + :attr:`StackLattice.central_axis`, :attr:`StackLattice.base_coordinate`, + :attr:`StackLattice.pitch`, :attr:`StackLattice.outer`, + and :attr:`StackLattice.universes` properties need to be set. + + Most methods for this class use a natural indexing scheme wherein elements + are assigned an index corresponding to their position relative to the + :math:`x`, :math:`y`, or :math:`z` bases, depending on the lattice + orientation, i.e., an index of + (0) in the lattice oriented to the :math:`z` basis gives the element whose + `z` coordinate is the smallest. However, note that when universes are + assigned to lattice elements using the + :attr:`StackLattice.universes` property, the array indices do not + correspond to natural indices. + + Parameters + ---------- + lattice_id : int, optional + Unique identifier for the lattice. If not specified, an identifier will + automatically be assigned. + name : str, optional + Name of the lattice. If not specified, the name is the empty string. + + Attributes + ---------- + id : int + Unique identifier for the lattice + name : str + Name of the lattice + orientation : {'x', 'y', 'z'} + Lattice central axis. Defaults to 'z'. + is_uniform : bool + Marks if the lattice is uniform or not. Defaults to True. + pitch : float or iterable of float + Pitch of the lattice in cm. If an iterable of float, then the ith entry + of the iterable is the width of the ith lattice element in cm. + outer : openmc.Universe + A universe to fill all space outside the lattice + universes : Iterable of openmc.Universe + A one-dimensional list/array of universes filling each element of the + lattice + central_axis : Iterable of float + The :math:`(y,z)`, :math:`(x,z)`, or :math:`(x,y)` coordinates of the + central axis of the lattice, depending on the lattice orientation + base_coordinate : float + The coordinate of the base layer of the lattice. Subsequent layers + are stacked towards the positive side of the orientation axis. + indices : list of tuple + A list of all possible lattice element indices. These + indices correspond to indices in the + :attr:`StackLattice.universes`property. + num_layers : int + An integer representing the number of lattice + layers along the orientation axis. + + """ + + def __init__(self, lattice_id=None, name=''): + super().__init__(lattice_id, name) + self.ndim = 1 + + # Initialize Lattice class attributes + self._central_axis = None + self._is_uniform = True + self.orientation = 'z' + + def __repr__(self): + string = 'StackLattice\n' + string += '{0: <16}{1}{2}\n'.format('\tID', '=\t', self._id) + string += '{0: <16}{1}{2}\n'.format('\tName', '=\t', self._name) + string += '{0: <16}{1}{2}\n'.format('\tOrientation', '=\t', + self._orientation) + string += '{0: <16}{1}{2}\n'.format('\t# Layers', + '=\t', self.num_layers) + string += '{0: <16}{1}{2}\n'.format('\tUniform', + '=\t', self._is_uniform) + string += '{0: <16}{1}{2}\n'.format('\tPitch', '=\t', self.pitch) + string += '{0: <16}{1}{2}\n'.format('\tbase_coordinate', '=\t', + self._base_coordinate) + string += '{0: <16}{1}{2}\n'.format('\tcentral_axis', '=\t', + self._central_axis) + if self._outer is not None: + string += '{0: <16}{1}{2}\n'.format('\tOuter', '=\t', + self._outer._id) + else: + string += '{0: <16}{1}{2}\n'.format('\tOuter', '=\t', + self._outer) + + string += '{: <16}\n'.format('\tUniverses, ' + + f'{self.orientation}-coord') + + for i, universe in enumerate(np.ravel(self._universes)): + string += f'{universe._id}, {self._layer_boundaries[i]}' + string += '\n' + + string = string.rstrip('\n') + + return string + + @property + def num_layers(self): + if self.universes is not None: + return len(self.universes) + else: + raise ValueError('Number of layers cannot be determined until ' + 'the universes has been set.') + + @property + def central_axis(self): + return self._central_axis + + @property + def base_coordinate(self): + return self._base_coordinate + + @property + def orientation(self): + return self._orientation + + @property + def is_uniform(self): + return self._is_uniform + + @property + def indices(self): + return list(np.ogrid[:self.num_layers]) + + @property + def _natural_indices(self): + """Iterate over all possible lattice element indices. + + """ + n = self.num_layers + for i in range(n): + yield i + + @central_axis.setter + def central_axis(self, central_axis): + cv.check_type('lattice central_axis', central_axis, Iterable, Real) + cv.check_length('lattice central_axis', central_axis, 2) + self._central_axis = central_axis + + @base_coordinate.setter + def base_coordinate(self, base_coordinate): + cv.check_type('lattice base_layer_coordinate', base_coordinate, Real) + self._base_coordinate = base_coordinate + + @orientation.setter + def orientation(self, orientation): + cv.check_value('orientation', orientation.lower(), ('x', 'y', 'z')) + + if orientation == 'x': + self._orientation_idx = 0 + elif orientation == 'y': + self._orientation_idx = 1 + else: + self._orientation_idx = 2 + self._orientation = orientation.lower() + + @is_uniform.setter + def is_uniform(self, is_uniform): + cv.check_type('is_uniform', is_uniform, bool) + self._is_uniform = is_uniform + + @Lattice.pitch.setter + def pitch(self, pitch): + if np.all(self.universes[0] == self.universes) and self._is_uniform: + cv.check_type('lattice pitch', pitch, Real) + _layer_boundaries = pitch * np.arange(1, self.num_layers + 1) + else: + cv.check_type('lattice pitch', pitch, Iterable, Real) + cv.check_length( + 'lattice pitch', + pitch, + self.num_layers, + self.num_layers) + _layer_boundaries = [pitch[0]] + for p in pitch[1:]: + _layer_boundaries += [_layer_boundaries[-1] + p] + _layer_boundaries = np.asarray(_layer_boundaries) + + _layer_boundaries += self._base_coordinate + self._layer_boundaries = np.insert( + _layer_boundaries, 0, self._base_coordinate) + self._pitch = pitch + + @Lattice.universes.setter + def universes(self, universes): + cv.check_iterable_type( + 'lattice universes', + universes, + openmc.UniverseBase, + min_depth=1, + max_depth=1) + universes = np.asarray(universes) + if len(np.shape(universes)) > 1: + raise SyntaxWarning( + "StackLattice universe array should be 1 dimensional.") + + self._universes = universes + cv.check_length('lattice universes', universes, self.num_layers) + + def find_element(self, point): + """Determine index of lattice element and local coordinates for a point + + Parameters + ---------- + point : Iterable of float + Cartesian coordinates of point + + Returns + ------- + int + The corresponding lattice element index + 3-tuple of float + Cartesian coordinates of the point in the corresponding lattice + element coordinate system + + """ + # find the level: + p = point[self._orientation_idx] + idx = 0 + if self._is_uniform: + idx = floor((p - self.base_coordinate) / self.pitch) + else: + if (p < self.base_coordinate): + idx = -1 # -1 index, not the last index + elif (p >= self._layer_boundaries[-1]): + idx = self.num_layers + else: + while not( + p >= self._layer_boundaries[idx] and p < self._layer_boundaries[idx + 1]): + idx += 1 + + return idx, self.get_local_coordinates(point, idx) + + def get_local_coordinates(self, point, idx): + """Determine local coordinates of a point within a lattice element + + Parameters + ---------- + point : Iterable of float + Cartesian coordinates of point + idx : int + index of lattice element + + Returns + ------- + 3-tuple of float + Cartesian coordinates of point in the lattice element coordinate + system + """ + x, y, z = point + c1, c2 = self.central_axis + + if self.orientation == 'x': + y -= c1 + z -= c2 + elif self.orientation == 'y': + x -= c1 + z -= c2 + else: + x -= c1 + y -= c2 + + local_point = [x, y, z] + if idx < 0: + local_point[self._orientation_idx] -= self._base_coordinate + elif idx >= self.num_layers: + local_point[self._orientation_idx] -= self._layer_boundaries[-2] + else: + local_point[self._orientation_idx] -= self._layer_boundaries[idx] + + return tuple(local_point) + + def get_universe_index(self, idx): + """Return index in the universes array corresponding + to a lattice element index + + Parameters + ---------- + idx : int + Lattice element index + + Returns + ------- + int + Index used when setting the :attr:`StackLattice.universes` property + + """ + return idx + + def is_valid_index(self, idx): + """Determine whether lattice element index is within defined range + + Parameters + ---------- + idx : int + Lattice element index + + Returns + ------- + bool + Whether index is valid + + """ + return (0 <= idx < self.num_layers) + + def create_xml_subelement(self, xml_element, memo=None): + """Add the lattice xml representation to an incoming xml element + + Parameters + ---------- + xml_element : xml.etree.ElementTree.Element + XML element to be added to + + memo : set or None + A set of object id's representing geometry entities already + written to the xml_element. This parameter is used internally + and should not be specified by users. + + Returns + ------- + None + + """ + if memo and self in memo: + return + if memo is not None: + memo.add(self) + + lattice_subelement = ET.Element("stack_lattice") + lattice_subelement.set("id", str(self._id)) + if len(self._name) > 0: + lattice_subelement.set("name", str(self._id)) + + # Export the Lattice cell pitch + pitch = ET.SubElement(lattice_subelement, "pitch") + if self._is_uniform: + pitch.text = str(self._pitch) + else: + pitch.text = ' '.join(map(str, self._pitch)) + + # Export Lattice Uniformity + uniform = ET.SubElement(lattice_subelement, "is_uniform") + uniform.text = str(self._is_uniform).lower() + + # Export the Lattice outer Universe (if specified) + if self._outer is not None: + outer = ET.SubElement(lattice_subelement, "outer") + outer.text = str(self._outer._id) + self._outer.create_xml_subelement(xml_element, memo) + + # Export numer of Lattice layers + num_layers = ET.SubElement(lattice_subelement, "num_layers") + num_layers.text = str(self.num_layers) + + # Export lattice orientation + lattice_subelement.set("orientation", self._orientation) + + # Export Lattice central axis + central_axis = ET.SubElement(lattice_subelement, "central_axis") + central_axis.text = ' '.join(map(str, self._central_axis)) + + # Export Lattice base coordinate + base_coordinate = ET.SubElement(lattice_subelement, "base_coordinate") + base_coordinate.text = str(self._base_coordinate) + + # Export the Lattice nested Universe IDs - column major for Fortran + universe_ids = '\n' + + # 1D stack + for l in range(self.num_layers): + universe = self._universes[l] + # Append Universe ID to the Lattice XML subelement + universe_ids += f'{universe._id} ' + + # Create XML subelement for this Universe + universe.create_xml_subelement(xml_element, memo) + + # Remove trailing newline character from Universe IDs string + universe_ids = universe_ids.rstrip('\n') + + universes = ET.SubElement(lattice_subelement, "universes") + universes.text = universe_ids + + # Append the XML subelement for this Lattice to the XML element + xml_element.append(lattice_subelement) + + @classmethod + def from_xml_element(cls, elem, get_universe): + """Generate nonuniform stack lattice from XML element + + Parameters + ---------- + elem : xml.etree.ElementTree.Element + `` element + get_universe : function + Function returning universe (defined in + :meth:`openmc.Geometry.from_xml`) + + Returns + ------- + StackLattice + Stack lattice + + """ + lat_id = int(get_text(elem, 'id')) + name = get_text(elem, 'name') + orientation = get_text(elem, 'orientation') + lat = cls(lat_id, name) + lat.orientation = orientation + lat.central_axis = [float(i) + for i in get_text(elem, 'central_axis').split()] + lat.base_coordinate = float(get_text(elem, 'base_coordinate')) + lat.is_uniform = bool(get_text(elem, 'is_uniform')) + lat.pitch = [float(i) for i in get_text(elem, 'pitch').split()] + if len(lat.pitch) == 1: + lat.pitch = lat.pitch[0] + outer = get_text(elem, 'outer') + if outer is not None: + lat.outer = get_universe(int(outer)) + + # Get array of universes + num_layers = get_text(elem, 'num_layers') + uarray = np.array([get_universe(int(i)) for i in + get_text(elem, 'universes').split()]) + uarray.shape = int(num_layers) + lat.universes = uarray + return lat + + @classmethod + def from_hdf5(cls, group, universes): + """Create nonuniform stack lattice from HDF5 group + + Parameters + ---------- + group : h5py.Group + Group in HDF5 file + universes : dict + Dictionary mapping universe IDs to instances of + :class:`openmc.Universe`. + + Returns + ------- + openmc.StackLattice + Stack lattice + + """ + + num_layers = group['num_layers'][()] + central_axis = group['central_axis'][...] + base_coordinate = group['base_coordinate'][()] + is_uniform = bool(group['is_uniform'][...]) + pitch = group['pitch'][...] + if len(pitch) == 1: + pitch = pitch[0] + outer = group['outer'][()] + universe_ids = group['universes'][...].flatten() + + if 'orientation' in group: + orientation = group['orientation'][()].decode() + else: + orientation = "z" + + # Create the Lattice + lattice_id = int(group.name.split('/')[-1].lstrip('lattice ')) + name = group['name'][()].decode() if 'name' in group else '' + lattice = cls(lattice_id, name) + lattice.central_axis = central_axis + lattice.base_coordinate = base_coordinate + lattice.orientation = orientation + lattice.is_uniform = is_uniform + + # If the Universe specified outer the Lattice is not void + if outer >= 0: + lattice.outer = universes[outer] + + # Build array of Universe pointers for the Lattice + uarray = np.empty(universe_ids.shape, dtype=openmc.Universe) + + for l in range(num_layers): + uarray[l] = universes[universe_ids[l]] + + # Set the universes for the lattice + lattice.universes = uarray + lattice.pitch = pitch + + return lattice diff --git a/src/geometry.cpp b/src/geometry.cpp index 216650ff189..fe55b2b4178 100644 --- a/src/geometry.cpp +++ b/src/geometry.cpp @@ -390,10 +390,11 @@ BoundaryInfo distance_to_boundary(Particle& p) // also means the lat.type attribute can be removed) std::pair> lattice_distance; switch (lat.type_) { - case LatticeType::rect: + case LatticeType::rect: { lattice_distance = lat.distance(r, u, coord.lattice_i); break; - case LatticeType::hex: + } + case LatticeType::hex: { auto& cell_above {model::cells[p.coord(i - 1).cell]}; Position r_hex {p.coord(i - 1).r}; r_hex -= cell_above->translation_; @@ -404,6 +405,12 @@ BoundaryInfo distance_to_boundary(Particle& p) lattice_distance = lat.distance(r_hex, u, coord.lattice_i); break; } + case LatticeType::stack: { + lattice_distance = lat.distance(r, u, coord.lattice_i); + break; + } + } + d_lat = lattice_distance.first; level_lat_trans = lattice_distance.second; diff --git a/src/hdf5_interface.cpp b/src/hdf5_interface.cpp index d53ad38790b..0eae17db6b5 100644 --- a/src/hdf5_interface.cpp +++ b/src/hdf5_interface.cpp @@ -645,6 +645,19 @@ void write_int(hid_t group_id, int ndim, const hsize_t* dims, const char* name, group_id, ndim, dims, name, H5T_NATIVE_INT, H5S_ALL, indep, buffer); } +void write_bool(hid_t group_id, int ndim, const hsize_t* dims, const char* name, + const hbool_t* buffer, bool indep) +{ + write_dataset_lowlevel( + group_id, ndim, dims, name, H5T_NATIVE_HBOOL, H5S_ALL, indep, buffer); +} + +void write_bool(hid_t group_id, const char* name, const bool b, bool indep) +{ + const hbool_t hb = (hbool_t)b; + write_bool(group_id, 0, nullptr, name, &hb, indep); +} + void write_llong(hid_t group_id, int ndim, const hsize_t* dims, const char* name, const long long* buffer, bool indep) { diff --git a/src/lattice.cpp b/src/lattice.cpp index fa2e2828ecf..44ccb56113e 100644 --- a/src/lattice.cpp +++ b/src/lattice.cpp @@ -4,12 +4,14 @@ #include #include +#include #include "openmc/cell.h" #include "openmc/error.h" #include "openmc/geometry.h" #include "openmc/geometry_aux.h" #include "openmc/hdf5_interface.h" +#include "openmc/settings.h" #include "openmc/string_utils.h" #include "openmc/vector.h" #include "openmc/xml_interface.h" @@ -514,7 +516,7 @@ HexLattice::HexLattice(pugi::xml_node lat_node) : Lattice {lat_node} if (univ_words.size() != n_univ) { fatal_error(fmt::format( "Expected {} universes for a hexagonal lattice with {} rings and {} " - "axial levels but {} were specified.", + "axial layers but {} were specified.", n_univ, n_rings_, n_axial_, univ_words.size())); } @@ -1078,6 +1080,347 @@ void HexLattice::to_hdf5_inner(hid_t lat_group) const write_int(lat_group, 3, dims, "universes", out.data(), false); } +//============================================================================== +// StackLattice implementation +//============================================================================== + +StackLattice::StackLattice(pugi::xml_node lat_node) : Lattice {lat_node} +{ + type_ = LatticeType::stack; + + // Read the lattice orientation. Default to 'z'. + if (check_for_node(lat_node, "orientation")) { + std::string orientation = get_node_value(lat_node, "orientation"); + if (orientation == "z") { + orientation_ = Orientation::z; + is_3d_ = true; // only in the sense that we intersect the z-plane + } else if (orientation == "y") { + orientation_ = Orientation::y; + } else if (orientation == "x") { + orientation_ = Orientation::x; + } else { + fatal_error("Unrecognized orientation '" + orientation + + "' for lattice " + std::to_string(id_)); + } + } else { + orientation_ = Orientation::z; + } + + // Read the number of lattice cells along its axis. + std::string n_layers_str {get_node_value(lat_node, "num_layers")}; + vector n_layers_words {split(n_layers_str)}; + if (n_layers_words.size() == 1) { + n_layers_ = std::stoi(n_layers_words[0]); + if (orientation_ == Orientation::z) { + orientation_idx_ = 2; + } else if (orientation_ == Orientation::y) { + orientation_idx_ = 1; + } else { + orientation_idx_ = 0; + } + } else { + fatal_error("Stack lattice must be one dimensional."); + } + + // Read the lattice central-axis location. + std::string ca_str {get_node_value(lat_node, "central_axis")}; + vector ca_words {split(ca_str)}; + if (ca_words.size() != 2) { + fatal_error("Number of entries on must be 2."); + } + central_axis_[0] = stod(ca_words[0]); + central_axis_[1] = stod(ca_words[1]); + + // Read the lattice base coordinate + std::string base_coord_str {get_node_value(lat_node, "base_coordinate")}; + vector base_coord_words {split(base_coord_str)}; + if (base_coord_words.size() != 1) { + fatal_error("Number of entries on must be 1."); + } + base_coordinate_ = stod(base_coord_words[0]); + + // Read if the lattice is uniform + std::string is_uniform_str {get_node_value(lat_node, "is_uniform")}; + if (is_uniform_str == "true") { + is_uniform_ = true; + } else { + is_uniform_ = false; + } + + // Read the lattice pitches. + std::string pitch_str {get_node_value(lat_node, "pitch")}; + vector pitch_words {split(pitch_str)}; + if (pitch_words.size() == 1 && is_uniform_) { + pitch_.push_back(stod(pitch_words[0])); + } else if (pitch_words.size() != n_layers_) { + fatal_error("Number of entries on must be the same as the " + "number of entries on ."); + } + + // Layers helper attribute + layer_boundaries_.push_back(base_coordinate_); + if (is_uniform_) { + for (int i = 1; i < n_layers_; i++) { + layer_boundaries_.push_back(layer_boundaries_[i - 1] + pitch_[0]); + } + layer_boundaries_.push_back(pitch_[0] + layer_boundaries_[n_layers_ - 1]); + + } else { + pitch_.push_back(stod(pitch_words[0])); + for (int i = 1; i < n_layers_; i++) { + pitch_.push_back(stod(pitch_words[i])); + layer_boundaries_.push_back(pitch_[i - 1] + layer_boundaries_[i - 1]); + } + layer_boundaries_.push_back( + pitch_[n_layers_ - 1] + layer_boundaries_[n_layers_ - 1]); + } + + // Read the universes and make sure the correct number was specified. + std::string univ_str {get_node_value(lat_node, "universes")}; + vector univ_words {split(univ_str)}; + if (univ_words.size() != n_layers_) { + fatal_error( + fmt::format("Expected {} universes for a stack lattice of size {} but {} " + "were specified.", + n_layers_, n_layers_, univ_words.size())); + } + + // Parse the universes. + universes_.resize(n_layers_, C_NONE); + for (int i = 0; i < n_layers_; i++) { + universes_[i] = std::stoi(univ_words[i]); + } +} + +//============================================================================== + +int32_t const& StackLattice::operator[](array const& i_xyz) +{ + return universes_[get_flat_index(i_xyz)]; +} + +//============================================================================== + +bool StackLattice::are_valid_indices(array const& i_xyz) const +{ + bool is_valid; + if (orientation_ == Orientation::x) { + is_valid = ((i_xyz[0] >= 0) && (i_xyz[0] < n_layers_) && (i_xyz[1] == 0) && + (i_xyz[2] == 0)); + } else if (orientation_ == Orientation::y) { + is_valid = ((i_xyz[0] == 0) && (i_xyz[1] >= 0) && (i_xyz[1] < n_layers_) && + (i_xyz[2] == 0)); + } else { + is_valid = ((i_xyz[0] == 0) && (i_xyz[1] == 0) && (i_xyz[2] >= 0) && + (i_xyz[2] < n_layers_)); + } + return is_valid; +} + +//============================================================================== + +std::pair> StackLattice::distance( + Position r, Direction u, const array& i_xyz) const +{ + + double c0, u_c, c; + array lattice_trans; + + // Determine the oncoming edge. + if (orientation_ == Orientation::x) { + u_c = u.x; + c = r.x; + lattice_trans = {1, 0, 0}; + } else if (orientation_ == Orientation::y) { + u_c = u.y; + c = r.y; + lattice_trans = {0, 1, 0}; + } else { + u_c = u.z; + c = r.z; + lattice_trans = {0, 0, 1}; + } + if (is_uniform_) { + c0 = copysign(pitch_[0], u_c); + } else { + c0 = copysign(pitch_[i_xyz[orientation_idx_]], u_c); + } + + // Top and bottom sides + double d {INFTY}; + if ((std::abs(c - c0) > FP_PRECISION) && u_c != 0) { + d = (c0 - c) / u_c; + if (d < 0) { + int i = 1; + } + if (u_c <= 0) { + lattice_trans[orientation_idx_] = -1; + } + } + + return {d, lattice_trans}; +} + +//============================================================================== + +void StackLattice::get_indices( + Position r, Direction u, array& result) const +{ + double r_c, u_c; + if (orientation_ == Orientation::x) { + r_c = r.x; + u_c = u.x; + result[1] = 0; + result[2] = 0; + } else if (orientation_ == Orientation::y) { + r_c = r.y; + u_c = u.y; + result[0] = 0; + result[2] = 0; + } else { + r_c = r.z; + u_c = u.z; + result[0] = 0; + result[1] = 0; + } + + // Determine axis index, accounting for coincidence + double ic_ = r_c - base_coordinate_; + + if (is_uniform_ || (r_c < base_coordinate_)) { + ic_ /= pitch_[0]; + } else { + // nonuniform case + if (r_c >= layer_boundaries_[n_layers_]) { + ic_ = n_layers_ + + ((r_c - layer_boundaries_[n_layers_]) / pitch_[n_layers_ - 1]); + } else { + ic_ = 0; + while (!((r_c >= layer_boundaries_[ic_]) and + (r_c < layer_boundaries_[ic_ + 1]))) { + ic_ += 1; + } + ic_ = (double)ic_ + ((r_c - layer_boundaries_[ic_]) / pitch_[ic_]); + } + } + if (settings::verbosity >= 10) { + write_message( + fmt::format( + " Lattice index {} in lattice {}. Current position r={}, u={}", ic_, + id_, r, u), + 1); + } + long ic_close {std::lround(ic_)}; + if (coincident(ic_, ic_close)) { + result[orientation_idx_] = (u_c > 0) ? ic_close : ic_close - 1; + } else { + result[orientation_idx_] = std::floor(ic_); + } +} + +int StackLattice::get_flat_index(const array& i_xyz) const +{ + return i_xyz[orientation_idx_]; +} + +//============================================================================== + +Position StackLattice::get_local_position( + Position r, const array& i_xyz) const +{ + double* c; + if (orientation_ == Orientation::x) { + r.y -= central_axis_[0]; + r.z -= central_axis_[1]; + c = &r.x; + } else if (orientation_ == Orientation::y) { + r.x -= central_axis_[0]; + r.z -= central_axis_[1]; + c = &r.y; + } else { + r.x -= central_axis_[0]; + r.y -= central_axis_[1]; + c = &r.z; + } + + if (i_xyz[orientation_idx_] < 0) { + *c -= layer_boundaries_[0]; + } else if (i_xyz[orientation_idx_] >= n_layers_) { + *c -= layer_boundaries_[n_layers_ - 2]; + } else { + *c -= layer_boundaries_[i_xyz[orientation_idx_]]; + } + return r; +} + +//============================================================================== + +int32_t& StackLattice::offset(int map, array const& i_xyz) +{ + return offsets_[n_layers_ * map + i_xyz[orientation_idx_]]; +} + +//============================================================================== + +int32_t StackLattice::offset(int map, int indx) const +{ + return offsets_[n_layers_ * map + indx]; +} + +//============================================================================== + +std::string StackLattice::index_to_string(int indx) const +{ + std::string out {std::to_string(indx)}; + return out; +} + +//============================================================================== + +void StackLattice::to_hdf5_inner(hid_t lat_group) const +{ + // Write basic lattice information. + write_string(lat_group, "type", "stack", false); + write_bool(lat_group, "is_uniform", is_uniform_, false); + write_dataset(lat_group, "pitch", pitch_); + write_dataset(lat_group, "central_axis", central_axis_); + write_dataset(lat_group, "base_coordinate", base_coordinate_); + write_dataset(lat_group, "num_layers", n_layers_); + + // Write lattice orientation + if (orientation_ == Orientation::x) { + write_string(lat_group, "orientation", "x", false); + } else if (orientation_ == Orientation::y) { + write_string(lat_group, "orientation", "y", false); + } else { + write_string(lat_group, "orientation", "z", false); + } + + // Write the universe ids. The convention here is that + // the lower-indexed universes are closer to the base coordinate + hsize_t nc {static_cast(n_layers_)}; + vector out(nc * 1 * 1); + hsize_t dims[3]; + for (int j = 0; j < nc; j++) { + out[j] = model::universes[universes_[j]]->id_; + } + if (orientation_ == Orientation::x) { + dims[0] = nc; + dims[1] = 1; + dims[2] = 1; + } else if (orientation_ == Orientation::y) { + dims[0] = 1; + dims[1] = nc; + dims[2] = 1; + } else { + dims[0] = 1; + dims[1] = 1; + dims[2] = nc; + } + + write_int(lat_group, 3, dims, "universes", out.data(), false); +} + //============================================================================== // Non-method functions //============================================================================== @@ -1090,6 +1433,9 @@ void read_lattices(pugi::xml_node node) for (pugi::xml_node lat_node : node.children("hex_lattice")) { model::lattices.push_back(make_unique(lat_node)); } + for (pugi::xml_node lat_node : node.children("stack_lattice")) { + model::lattices.push_back(make_unique(lat_node)); + } // Fill the lattice map. for (int i_lat = 0; i_lat < model::lattices.size(); i_lat++) { diff --git a/tests/regression_tests/lattice_stack_nonuniform/__init__.py b/tests/regression_tests/lattice_stack_nonuniform/__init__.py new file mode 100644 index 00000000000..e69de29bb2d diff --git a/tests/regression_tests/lattice_stack_nonuniform/inputs_true.dat b/tests/regression_tests/lattice_stack_nonuniform/inputs_true.dat new file mode 100644 index 00000000000..a9100a7f1ff --- /dev/null +++ b/tests/regression_tests/lattice_stack_nonuniform/inputs_true.dat @@ -0,0 +1,51 @@ + + + + + + + + + + + 1.0 2.0 1.0 2.0 1.0 2.0 1.0 2.0 1.0 2.0 1.0 2.0 1.0 2.0 1.0 2.0 1.0 2.0 1.0 2.0 1.0 2.0 1.0 2.0 1.0 2.0 1.0 2.0 1.0 2.0 1.0 2.0 1.0 2.0 1.0 2.0 1.0 2.0 1.0 2.0 1.0 2.0 1.0 2.0 1.0 2.0 1.0 2.0 1.0 2.0 1.0 2.0 1.0 2.0 1.0 2.0 1.0 2.0 1.0 2.0 1.0 2.0 1.0 2.0 1.0 2.0 1.0 2.0 1.0 2.0 1.0 2.0 1.0 2.0 1.0 2.0 1.0 2.0 1.0 2.0 1.0 2.0 1.0 2.0 1.0 2.0 1.0 2.0 1.0 2.0 1.0 2.0 1.0 2.0 1.0 2.0 1.0 2.0 1.0 2.0 1.0 2.0 1.0 2.0 1.0 2.0 1.0 2.0 1.0 2.0 1.0 2.0 1.0 2.0 1.0 2.0 1.0 2.0 1.0 2.0 1.0 2.0 1.0 2.0 1.0 2.0 1.0 2.0 1.0 2.0 1.0 2.0 1.0 2.0 1.0 2.0 1.0 2.0 1.0 2.0 1.0 2.0 1.0 2.0 1.0 2.0 1.0 2.0 1.0 2.0 1.0 2.0 1.0 2.0 1.0 2.0 1.0 2.0 1.0 2.0 1.0 2.0 1.0 2.0 1.0 2.0 1.0 2.0 1.0 2.0 1.0 2.0 1.0 2.0 1.0 2.0 1.0 2.0 1.0 2.0 1.0 2.0 1.0 2.0 1.0 2.0 1.0 2.0 1.0 2.0 1.0 2.0 1.0 2.0 1.0 2.0 1.0 2.0 1.0 2.0 + false + 200 + 0.0 0.0 + 0.0 + +1 2 1 2 1 2 1 2 1 2 1 2 1 2 1 2 1 2 1 2 1 2 1 2 1 2 1 2 1 2 1 2 1 2 1 2 1 2 1 2 1 2 1 2 1 2 1 2 1 2 1 2 1 2 1 2 1 2 1 2 1 2 1 2 1 2 1 2 1 2 1 2 1 2 1 2 1 2 1 2 1 2 1 2 1 2 1 2 1 2 1 2 1 2 1 2 1 2 1 2 1 2 1 2 1 2 1 2 1 2 1 2 1 2 1 2 1 2 1 2 1 2 1 2 1 2 1 2 1 2 1 2 1 2 1 2 1 2 1 2 1 2 1 2 1 2 1 2 1 2 1 2 1 2 1 2 1 2 1 2 1 2 1 2 1 2 1 2 1 2 1 2 1 2 1 2 1 2 1 2 1 2 1 2 1 2 1 2 1 2 1 2 1 2 1 2 1 2 1 2 + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + eigenvalue + 1000 + 10 + 5 + diff --git a/tests/regression_tests/lattice_stack_nonuniform/results_true.dat b/tests/regression_tests/lattice_stack_nonuniform/results_true.dat new file mode 100644 index 00000000000..c3f32c6e4fa --- /dev/null +++ b/tests/regression_tests/lattice_stack_nonuniform/results_true.dat @@ -0,0 +1,2 @@ +k-combined: +1.850801E+00 5.892615E-03 diff --git a/tests/regression_tests/lattice_stack_nonuniform/test.py b/tests/regression_tests/lattice_stack_nonuniform/test.py new file mode 100644 index 00000000000..41cdc31fa75 --- /dev/null +++ b/tests/regression_tests/lattice_stack_nonuniform/test.py @@ -0,0 +1,77 @@ +from tests.testing_harness import PyAPITestHarness + +import numpy as np +import openmc +import pytest + + +def nonuniform_stack_lattice_model(): + model = openmc.model.Model() + + uo2 = openmc.Material(name='UO2') + uo2.set_density('g/cm3', 10.0) + uo2.add_nuclide('U235', 1.0) + uo2.add_nuclide('O16', 2.0) + water = openmc.Material(name='light water') + water.add_nuclide('H1', 2.0) + water.add_nuclide('O16', 1.0) + water.set_density('g/cm3', 1.0) + water.add_s_alpha_beta('c_H_in_H2O') + model.materials.extend([uo2, water]) + + rc = 0.4 + h1 = 1.0 + h2 = 2.0 + h_avg = (h1 + h2) / 2 + cyl = openmc.ZCylinder(r=rc) + top1 = openmc.ZPlane(z0=h1) + top2 = openmc.ZPlane(z0=h2) + bottom = openmc.ZPlane(z0=0.) + pellet1 = -cyl & -top1 & +bottom + pellet2 = -cyl & -top2 & +bottom + water_slice1 = +cyl & -top1 & +bottom + water_slice2 = +cyl & -top2 & +bottom + + fuel1 = openmc.Cell(fill=uo2, region=pellet1) + fuel2 = openmc.Cell(fill=uo2, region=pellet2) + water_reflector1 = openmc.Cell(fill=water, region=water_slice1) + water_reflector2 = openmc.Cell(fill=water, region=water_slice2) + layer1 = openmc.Universe(cells=[fuel1, water_reflector1]) + layer2 = openmc.Universe(cells=[fuel2, water_reflector2]) + + n_pellets = 200 + + top = openmc.ZPlane(z0=n_pellets * h_avg) + tb_refl = openmc.Cell(fill=water, region=-bottom | +top) + + d = 1.5 * rc + box = openmc.model.RectangularParallelepiped(-d, d, -d, d, 0. - d, + n_pellets * h_avg + d, + boundary_type='reflective') + + univs = [layer1, layer2] * int(n_pellets / 2) + pellet_stack = openmc.StackLattice() + pellet_stack.central_axis = (0., 0.) + pellet_stack.base_coordinate = 0. + pellet_stack.universes = univs + pellet_stack.is_uniform = False + pellet_stack.pitch = [h1, h2] * int(n_pellets / 2) + + stack_cell = openmc.Cell(fill=pellet_stack) + + pin_univ = openmc.Universe(cells=[stack_cell, tb_refl]) + + main_cell = openmc.Cell(fill=pin_univ, region=-box) + model.geometry = openmc.Geometry([main_cell]) + + model.settings.batches = 10 + model.settings.inactive = 5 + model.settings.particles = 1000 + + return model + + +def test_lattice_stack_nonuniform(): + model = nonuniform_stack_lattice_model() + harness = PyAPITestHarness('statepoint.10.h5', model) + harness.main() diff --git a/tests/regression_tests/lattice_stack_uniform/__init__.py b/tests/regression_tests/lattice_stack_uniform/__init__.py new file mode 100644 index 00000000000..e69de29bb2d diff --git a/tests/regression_tests/lattice_stack_uniform/inputs_true.dat b/tests/regression_tests/lattice_stack_uniform/inputs_true.dat new file mode 100644 index 00000000000..425b087e64f --- /dev/null +++ b/tests/regression_tests/lattice_stack_uniform/inputs_true.dat @@ -0,0 +1,48 @@ + + + + + + + + + 1.5 + true + 200 + 0.0 0.0 + 0.0 + +1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 + + + + + + + + + + + + + + + + + + + + + + + + + + + + + eigenvalue + 1000 + 10 + 5 + diff --git a/tests/regression_tests/lattice_stack_uniform/results_true.dat b/tests/regression_tests/lattice_stack_uniform/results_true.dat new file mode 100644 index 00000000000..068459959c3 --- /dev/null +++ b/tests/regression_tests/lattice_stack_uniform/results_true.dat @@ -0,0 +1,2 @@ +k-combined: +1.839870E+00 6.141508E-03 diff --git a/tests/regression_tests/lattice_stack_uniform/test.py b/tests/regression_tests/lattice_stack_uniform/test.py new file mode 100644 index 00000000000..3ec5155f300 --- /dev/null +++ b/tests/regression_tests/lattice_stack_uniform/test.py @@ -0,0 +1,69 @@ +import numpy as np +import openmc +import pytest + +from tests.testing_harness import PyAPITestHarness + + +def uniform_stack_lattice_model(): + model = openmc.model.Model() + + uo2 = openmc.Material(name='UO2') + uo2.set_density('g/cm3', 10.0) + uo2.add_nuclide('U235', 1.0) + uo2.add_nuclide('O16', 2.0) + water = openmc.Material(name='light water') + water.add_nuclide('H1', 2.0) + water.add_nuclide('O16', 1.0) + water.set_density('g/cm3', 1.0) + water.add_s_alpha_beta('c_H_in_H2O') + model.materials.extend([uo2, water]) + + rc = 0.4 + h = 1.5 + cyl = openmc.ZCylinder(r=rc) + top = openmc.ZPlane(z0=h) + bottom = openmc.ZPlane(z0=0.) + pellet = -cyl & -top & +bottom + water_slice = +cyl & -top & +bottom + + fuel = openmc.Cell(fill=uo2, region=pellet) + water_reflector = openmc.Cell(fill=water, region=water_slice) + layer = openmc.Universe(cells=[fuel, water_reflector]) + + n_pellets = 200 + + top = openmc.ZPlane(z0=n_pellets * h) + tb_refl = openmc.Cell(fill=water, region=-bottom | +top) + + univs = [layer] * n_pellets + pellet_stack = openmc.StackLattice() + pellet_stack.central_axis = (0., 0.) + pellet_stack.base_coordinate = 0. + pellet_stack.universes = univs + pellet_stack.is_uniform = True + pellet_stack.pitch = h + + stack_cell = openmc.Cell(fill=pellet_stack) + + pin_univ = openmc.Universe(cells=[stack_cell, tb_refl]) + + d = 1.5 * rc + box = openmc.model.RectangularParallelepiped(-d, d, -d, d, + 0. - d, n_pellets * h + d, + boundary_type='reflective') + + main_cell = openmc.Cell(fill=pin_univ, region=-box) + model.geometry = openmc.Geometry([main_cell]) + + model.settings.batches = 10 + model.settings.inactive = 5 + model.settings.particles = 1000 + + return model + + +def test_lattice_stack_uniform(): + model = uniform_stack_lattice_model() + harness = PyAPITestHarness('statepoint.10.h5', model) + harness.main() diff --git a/tests/unit_tests/test_lattice.py b/tests/unit_tests/test_lattice.py index 99bc284b0eb..b13f8c57131 100644 --- a/tests/unit_tests/test_lattice.py +++ b/tests/unit_tests/test_lattice.py @@ -45,7 +45,7 @@ def rlat2(pincell1, pincell2, uo2, water, zr): n = 3 u1, u2 = pincell1, pincell2 lattice = openmc.RectLattice() - lattice.lower_left = (-pitch*n/2, -pitch*n/2) + lattice.lower_left = (-pitch * n / 2, -pitch * n / 2) lattice.pitch = (pitch, pitch) lattice.outer = openmc.Universe(cells=[all_zr]) lattice.universes = [ @@ -77,7 +77,7 @@ def rlat3(pincell1, pincell2, uo2, water, zr): n = 3 u1, u2 = pincell1, pincell2 lattice = openmc.RectLattice() - lattice.lower_left = (-pitch*n/2, -pitch*n/2, -10.0) + lattice.lower_left = (-pitch * n / 2, -pitch * n / 2, -10.0) lattice.pitch = (pitch, pitch, 10.0) lattice.outer = openmc.Universe(cells=[all_zr]) lattice.universes = [ @@ -155,7 +155,53 @@ def hlat3(pincell1, pincell2, uo2, water, zr): return lattice -def test_get_nuclides(rlat2, rlat3, hlat2, hlat3): +@pytest.fixture(scope='module') +def slatU(pincell1, pincell2, uo2, water, zr): + """Uniform Stack lattice for testing.""" + + all_zr = openmc.Cell(fill=zr) + pitch = 1.2 + u1 = pincell1 + lattice = openmc.StackLattice() + lattice.orientation = 'z' + lattice.central_axis = (0., 0.) + lattice.base_coordinate = 0. + lattice.is_uniform = True + lattice.universes = [u1, u1, u1, u1] + lattice.pitch = pitch + lattice.outer = openmc.Universe(cells=[all_zr]) + + # Add extra attributes for comparison purpose + lattice.cells = [u1.fuel, u1.moderator, all_zr] + lattice.mats = [uo2, water, zr] + lattice.univs = [u1, lattice.outer] + return lattice + + +@pytest.fixture(scope='module') +def slatNU(pincell1, pincell2, uo2, water, zr): + """Nonuniform Stack lattice for testing.""" + + all_zr = openmc.Cell(fill=zr) + pitch = 1.2 + u1, u2 = pincell1, pincell2 + lattice = openmc.StackLattice() + lattice.orientation = 'z' + lattice.central_axis = (0., 0.) + lattice.base_coordinate = 0. + lattice.is_uniform = False + lattice.universes = [u1, u2, u1, u2] + lattice.pitch = [pitch, 1.5 * pitch, pitch, 1.5 * pitch] + lattice.outer = openmc.Universe(cells=[all_zr]) + + # Add extra attributes for comparison purpose + lattice.cells = [u1.fuel, u1.moderator, u2.fuel, u2.moderator, all_zr] + lattice.mats = [uo2, water, zr] + lattice.univs = [u1, u2, lattice.outer] + return lattice + + +def test_get_nuclides(rlat2, rlat3, hlat2, hlat3, slatU, slatNU): for lat in (rlat2, hlat2): nucs = lat.get_nuclides() assert sorted(nucs) == ['H1', 'O16', 'U235', @@ -164,27 +210,31 @@ def test_get_nuclides(rlat2, rlat3, hlat2, hlat3): nucs = lat.get_nuclides() assert sorted(nucs) == ['H1', 'H2', 'O16', 'U235', 'Zr90', 'Zr91', 'Zr92', 'Zr94', 'Zr96'] + for lat in (slatU, slatNU): + nucs = lat.get_nuclides() + assert sorted(nucs) == ['H1', 'O16', 'U235', + 'Zr90', 'Zr91', 'Zr92', 'Zr94', 'Zr96'] -def test_get_all_cells(rlat2, rlat3, hlat2, hlat3): - for lat in (rlat2, rlat3, hlat2, hlat3): +def test_get_all_cells(rlat2, rlat3, hlat2, hlat3, slatU, slatNU): + for lat in (rlat2, rlat3, hlat2, hlat3, slatU, slatNU): cells = set(lat.get_all_cells().values()) assert not cells ^ set(lat.cells) -def test_get_all_materials(rlat2, rlat3, hlat2, hlat3): - for lat in (rlat2, rlat3, hlat2, hlat3): +def test_get_all_materials(rlat2, rlat3, hlat2, hlat3, slatU, slatNU): + for lat in (rlat2, rlat3, hlat2, hlat3, slatU, slatNU): mats = set(lat.get_all_materials().values()) assert not mats ^ set(lat.mats) -def test_get_all_universes(rlat2, rlat3, hlat2, hlat3): - for lat in (rlat2, rlat3, hlat2, hlat3): +def test_get_all_universes(rlat2, rlat3, hlat2, hlat3, slatU, slatNU): + for lat in (rlat2, rlat3, hlat2, hlat3, slatU, slatNU): univs = set(lat.get_all_universes().values()) assert not univs ^ set(lat.univs) -def test_get_universe(rlat2, rlat3, hlat2, hlat3): +def test_get_universe(rlat2, rlat3, hlat2, hlat3, slatU, slatNU): u1, u2, outer = rlat2.univs assert rlat2.get_universe((0, 0)) == u2 assert rlat2.get_universe((1, 0)) == u1 @@ -217,8 +267,20 @@ def test_get_universe(rlat2, rlat3, hlat2, hlat3): assert hlat3.get_universe((0, -2, 0)) == u1 assert hlat3.get_universe((0, -2, 1)) == u3 + u1, outer = slatU.univs + assert slatU.get_universe(0) == u1 + assert slatU.get_universe(1) == u1 + assert slatU.get_universe(2) == u1 + assert slatU.get_universe(3) == u1 + + u1, u2, outer = slatNU.univs + assert slatNU.get_universe(0) == u1 + assert slatNU.get_universe(1) == u2 + assert slatNU.get_universe(2) == u1 + assert slatNU.get_universe(3) == u2 + -def test_find(rlat2, rlat3, hlat2, hlat3): +def test_find(rlat2, rlat3, hlat2, hlat3, slatU, slatNU): pitch = rlat2.pitch[0] seq = rlat2.find((0., 0., 0.)) assert seq[-1] == rlat2.cells[0] @@ -226,7 +288,7 @@ def test_find(rlat2, rlat3, hlat2, hlat3): assert seq[-1] == rlat2.cells[2] seq = rlat2.find((0., -pitch, 0.)) assert seq[-1] == rlat2.cells[0] - seq = rlat2.find((pitch*100, 0., 0.)) + seq = rlat2.find((pitch * 100, 0., 0.)) assert seq[-1] == rlat2.cells[-1] seq = rlat3.find((-pitch, pitch, 5.0)) assert seq[-1] == rlat3.cells[-2] @@ -236,7 +298,7 @@ def test_find(rlat2, rlat3, hlat2, hlat3): assert seq[-1] == hlat2.cells[2] seq = hlat2.find((0.5, 0., 0.)) assert seq[-1] == hlat2.cells[3] - seq = hlat2.find((sqrt(3)*pitch, 0., 0.)) + seq = hlat2.find((sqrt(3) * pitch, 0., 0.)) assert seq[-1] == hlat2.cells[0] seq = hlat2.find((0., pitch, 0.)) assert seq[-1] == hlat2.cells[2] @@ -248,7 +310,7 @@ def test_find(rlat2, rlat3, hlat2, hlat3): assert seq[-1] == hlat3.cells[2] seq = hlat3.find((0., -pitch, -5.)) assert seq[-1] == hlat3.cells[0] - seq = hlat3.find((sqrt(3)*pitch, 0., -5.)) + seq = hlat3.find((sqrt(3) * pitch, 0., -5.)) assert seq[-1] == hlat3.cells[0] # top of 3D lattice @@ -258,16 +320,47 @@ def test_find(rlat2, rlat3, hlat2, hlat3): assert seq[-1] == hlat3.cells[0] seq = hlat3.find((0., -pitch, 5.)) assert seq[-1] == hlat3.cells[-2] - seq = hlat3.find((sqrt(3)*pitch, 0., 5.)) + seq = hlat3.find((sqrt(3) * pitch, 0., 5.)) assert seq[-1] == hlat3.cells[0] - -def test_clone(rlat2, hlat2, hlat3): + pitch = slatU.pitch + seq = slatU.find((0., 0., 0.)) + assert seq[-1] == slatU.cells[0] + seq = slatU.find((0., 0., pitch)) + assert seq[-1] == slatU.cells[0] + seq = slatU.find((0., 0., 2 * pitch)) + assert seq[-1] == slatU.cells[0] + seq = slatU.find((0., 0., 3 * pitch)) + assert seq[-1] == slatU.cells[0] + + seq = slatNU.find((0., 0., 0.)) + assert seq[-1] == slatNU.cells[0] + seq = slatNU.find((0., 0., 2.9)) + assert seq[-1] == slatNU.cells[2] + seq = slatNU.find((0., 0., 4.1)) + assert seq[-1] == slatNU.cells[0] + seq = slatNU.find((0., 0., 5.9)) + assert seq[-1] == slatNU.cells[2] + + +def test_clone(rlat2, hlat2, hlat3, slatU, slatNU): rlat_clone = rlat2.clone() assert rlat_clone.id != rlat2.id assert rlat_clone.lower_left == rlat2.lower_left assert rlat_clone.pitch == rlat2.pitch + slatU_clone = slatU.clone() + assert slatU_clone.id != slatU.id + assert slatU_clone.central_axis == slatU.central_axis + assert slatU_clone.base_coordinate == slatU.base_coordinate + assert slatU_clone.pitch == slatU.pitch + + slatNU_clone = slatNU.clone() + assert slatU_clone.id != slatU.id + assert slatU_clone.central_axis == slatU.central_axis + assert slatU_clone.base_coordinate == slatU.base_coordinate + assert slatU_clone.pitch == slatU.pitch + hlat_clone = hlat2.clone() assert hlat_clone.id != hlat2.id assert hlat_clone.center == hlat2.center @@ -289,11 +382,13 @@ def test_clone(rlat2, hlat2, hlat3): assert c1.region == c2.region -def test_repr(rlat2, rlat3, hlat2, hlat3): +def test_repr(rlat2, rlat3, hlat2, hlat3, slatU, slatNU): repr(rlat2) repr(rlat3) repr(hlat2) repr(hlat3) + repr(slatU) + repr(slatNU) def test_indices_rect(rlat2, rlat3): @@ -331,6 +426,21 @@ def test_indices_hex(hlat2, hlat3): ) +def test_indices_stack(slatU, slatNU): + assert slatU.indices == [0, 1, 2, 3] + assert slatNU.indices == [0, 1, 2, 3] + slatU.orientation = 'x' + slatNU.orientation = 'x' + assert slatU.indices == [0, 1, 2, 3] + assert slatNU.indices == [0, 1, 2, 3] + slatU.orientation = 'y' + slatNU.orientation = 'y' + assert slatU.indices == [0, 1, 2, 3] + assert slatNU.indices == [0, 1, 2, 3] + slatU.orientation = 'z' + slatNU.orientation = 'z' + + def test_xml_rect(rlat2, rlat3): for lat in (rlat2, rlat3): geom = ET.Element('geometry') @@ -355,12 +465,37 @@ def test_xml_hex(hlat2, hlat3): assert len(elem.find('universes').text.split()) == len(lat.indices) +def test_xml_stack(slatU, slatNU): + for lat in (slatU, slatNU): + geom = ET.Element('geometry') + lat.create_xml_subelement(geom) + elem = geom.find('stack_lattice') + assert elem.tag == 'stack_lattice' + assert elem.get('id') == str(lat.id) + if lat.is_uniform: + assert len(elem.find('pitch').text.split()) == 1 + else: + assert len(elem.find('pitch').text.split()) == lat.num_layers + assert len(elem.find('universes').text.split()) == len(lat.indices) + + def test_show_indices(): for i in range(1, 11): lines = openmc.HexLattice.show_indices(i).split('\n') - assert len(lines) == 4*i - 3 + assert len(lines) == 4 * i - 3 lines_x = openmc.HexLattice.show_indices(i, 'x').split('\n') - assert len(lines_x) == 4*i - 3 + assert len(lines_x) == 4 * i - 3 + + +def test_unset_pitch(): + elem = ET.Element("dummy") + + stack_lattice = openmc.StackLattice() + stack_lattice.central_axis = (0., 0.) + stack_lattice.base_coordinate = 0. + + with pytest.raises(ValueError): + stack_lattice.create_xml_subelement(elem) def test_unset_universes():