crp/OpenMC
1
0
Fork 0
mirror of https://github.com/openmc-dev/openmc.git synced 2026-07-27 13:45:36 -04:00
Code Issues Projects Releases Packages Wiki Activity Actions

Merge branch 'develop' into cylinder-sector

Browse source
This commit is contained in:
yardasol 2022-04-07 18:55:51 -05:00 • committed by GitHub
commit 294f57318d
No known key found for this signature in database
GPG key ID: 4AEE18F83AFDEB23
20 changed files with 621 additions and 403 deletions
Download patch file Download diff file Expand all files Collapse all files

19
CMakeLists.txt
View file

@ -1,4 +1,4 @@
cmake_minimum_required(VERSION 3.3 FATAL_ERROR)
cmake_minimum_required(VERSION 3.10 FATAL_ERROR)
project(openmc C CXX)
# Set version numbers
@ -30,6 +30,7 @@ option(OPENMC_ENABLE_PROFILE "Compile with profiling flags"
option(OPENMC_ENABLE_COVERAGE "Compile with coverage analysis flags" OFF)
option(OPENMC_USE_DAGMC "Enable support for DAGMC (CAD) geometry" OFF)
option(OPENMC_USE_LIBMESH "Enable support for libMesh unstructured mesh tallies" OFF)
option(OPENMC_USE_MPI "Enable MPI" OFF)
#===============================================================================
# Set a default build configuration if not explicitly specified
@ -44,10 +45,8 @@ endif()
# MPI for distributed-memory parallelism
#===============================================================================
set(MPI_ENABLED FALSE)
if(${CMAKE_CXX_COMPILER} MATCHES "(mpi[^/]*|CC)$")
message(STATUS "Detected MPI wrapper: ${CMAKE_CXX_COMPILER}")
set(MPI_ENABLED TRUE)
if(OPENMC_USE_MPI)
find_package(MPI REQUIRED)
endif()
#===============================================================================
@ -108,7 +107,7 @@ endif()
find_package(HDF5 REQUIRED COMPONENTS C HL)
if(HDF5_IS_PARALLEL)
if(NOT MPI_ENABLED)
if(NOT OPENMC_USE_MPI)
message(FATAL_ERROR "Parallel HDF5 was detected, but the detected compiler,\
${CMAKE_CXX_COMPILER}, does not support MPI. An MPI-capable compiler must \
be used with parallel HDF5.")
@ -144,6 +143,7 @@ set(CMAKE_POSITION_INDEPENDENT_CODE ON)
if(OPENMC_ENABLE_PROFILE)
list(APPEND cxxflags -g -fno-omit-frame-pointer)
endif()
if(OPENMC_ENABLE_COVERAGE)
list(APPEND cxxflags --coverage)
list(APPEND ldflags --coverage)
@ -351,6 +351,7 @@ list(APPEND libopenmc_SOURCES
src/timer.cpp
src/thermal.cpp
src/track_output.cpp
src/universe.cpp
src/urr.cpp
src/volume_calc.cpp
src/weight_windows.cpp
@ -398,7 +399,7 @@ target_include_directories(libopenmc PRIVATE ${CMAKE_BINARY_DIR}/include)
if (HDF5_IS_PARALLEL)
target_compile_definitions(libopenmc PRIVATE -DPHDF5)
endif()
if (MPI_ENABLED)
if (OPENMC_USE_MPI)
target_compile_definitions(libopenmc PUBLIC -DOPENMC_MPI)
endif()
@ -440,6 +441,10 @@ if (PNG_FOUND)
target_link_libraries(libopenmc PNG::PNG)
endif()
if (OPENMC_USE_MPI)
target_link_libraries(libopenmc MPI::MPI_CXX)
endif()
#===============================================================================
# openmc executable
#===============================================================================

7
Dockerfile
View file

@ -59,8 +59,7 @@ ENV LIBMESH_INSTALL_DIR=$HOME/LIBMESH
ENV NJOY_REPO='https://github.com/njoy/NJOY2016'
# Setup environment variables for Docker image
ENV CC=/usr/bin/mpicc CXX=/usr/bin/mpicxx \
LD_LIBRARY_PATH=${DAGMC_INSTALL_DIR}/lib:$LD_LIBRARY_PATH \
ENV LD_LIBRARY_PATH=${DAGMC_INSTALL_DIR}/lib:$LD_LIBRARY_PATH \
OPENMC_CROSS_SECTIONS=/root/nndc_hdf5/cross_sections.xml \
OPENMC_ENDF_DATA=/root/endf-b-vii.1 \
DEBIAN_FRONTEND=noninteractive
@ -179,6 +178,7 @@ RUN mkdir -p ${HOME}/OpenMC && cd ${HOME}/OpenMC \
&& mkdir build && cd build ; \
if [ ${build_dagmc} = "on" ] && [ ${build_libmesh} = "on" ]; then \
cmake ../openmc \
-DOPENMC_USE_MPI=on \
-DHDF5_PREFER_PARALLEL=on \
-DOPENMC_USE_DAGMC=on \
-DOPENMC_USE_LIBMESH=on \
@ -186,18 +186,21 @@ RUN mkdir -p ${HOME}/OpenMC && cd ${HOME}/OpenMC \
fi ; \
if [ ${build_dagmc} = "on" ] && [ ${build_libmesh} = "off" ]; then \
cmake ../openmc \
-DOPENMC_USE_MPI=on \
-DHDF5_PREFER_PARALLEL=on \
-DOPENMC_USE_DAGMC=ON \
-DCMAKE_PREFIX_PATH=${DAGMC_INSTALL_DIR} ; \
fi ; \
if [ ${build_dagmc} = "off" ] && [ ${build_libmesh} = "on" ]; then \
cmake ../openmc \
-DOPENMC_USE_MPI=on \
-DHDF5_PREFER_PARALLEL=on \
-DOPENMC_USE_LIBMESH=on \
-DCMAKE_PREFIX_PATH=${LIBMESH_INSTALL_DIR} ; \
fi ; \
if [ ${build_dagmc} = "off" ] && [ ${build_libmesh} = "off" ]; then \
cmake ../openmc \
-DOPENMC_USE_MPI=on \
-DHDF5_PREFER_PARALLEL=on ; \
fi ; \
make 2>/dev/null -j${compile_cores} install \

4
cmake/OpenMCConfig.cmake.in
View file

@ -21,3 +21,7 @@ find_package(PNG)
if(NOT TARGET OpenMC::libopenmc)
include("${OpenMC_CMAKE_DIR}/OpenMCTargets.cmake")
endif()
if(@OPENMC_USE_MPI@)
find_package(MPI REQUIRED)
endif()

1
docs/source/pythonapi/model.rst
View file

@ -26,6 +26,7 @@ Composite Surfaces
:template: myclass.rst
openmc.model.CylinderSector
openmc.model.IsogonalOctagon
openmc.model.RectangularParallelepiped
openmc.model.RightCircularCylinder
openmc.model.XConeOneSided

26
docs/source/usersguide/install.rst
View file

@ -203,7 +203,7 @@ Prerequisites
respectively. To link against a parallel HDF5 library, make sure to set
the HDF5_PREFER_PARALLEL CMake option, e.g.::
CXX=mpicxx.mpich cmake -DHDF5_PREFER_PARALLEL=on ..
cmake -DHDF5_PREFER_PARALLEL=on -DOPENMC_USE_MPI=on ..
Note that the exact package names may vary depending on your particular
distribution and version.
@ -263,7 +263,7 @@ Prerequisites
installation should be specified as part of the ``CMAKE_PREFIX_PATH``
variable.::
CXX=mpicxx cmake -DOPENMC_USE_LIBMESH=on -DCMAKE_PREFIX_PATH=/path/to/libmesh/installation
cmake -DOPENMC_USE_LIBMESH=on -DOPENMC_USE_MPI=on -DCMAKE_PREFIX_PATH=/path/to/libmesh/installation
Note that libMesh is most commonly compiled with MPI support. If that
is the case, then OpenMC should be compiled with MPI support as well.
@ -352,6 +352,10 @@ OPENMC_ENABLE_COVERAGE
Compile and link code instrumented for coverage analysis. This is typically
used in conjunction with gcov_. (Default: off)
OPENMC_USE_MPI
Turns on compiling with MPI (default: off). For further information on MPI options,
please see the `FindMPI.cmake documentation <https://cmake.org/cmake/help/latest/module/FindMPI.html>`_.
To set any of these options (e.g., turning on profiling), the following form
should be used:
@ -384,24 +388,6 @@ Example of configuring for Debug mode:
cmake -DCMAKE_BUILD_TYPE=Debug /path/to/openmc
Compiling with MPI
++++++++++++++++++
To compile with MPI, set the :envvar:`CXX` environment variable to the path to
the MPI C++ wrapper. For example, in a bash shell:
.. code-block:: sh
export CXX=mpicxx
cmake /path/to/openmc
Note that in many shells, environment variables can be set for a single command,
i.e.
.. code-block:: sh
CXX=mpicxx cmake /path/to/openmc
Selecting HDF5 Installation
+++++++++++++++++++++++++++

62
include/openmc/cell.h
View file

@ -17,6 +17,7 @@
#include "openmc/neighbor_list.h"
#include "openmc/position.h"
#include "openmc/surface.h"
#include "openmc/universe.h"
#include "openmc/vector.h"
namespace openmc {
@ -48,36 +49,8 @@ namespace model {
extern std::unordered_map<int32_t, int32_t> cell_map;
extern vector<unique_ptr<Cell>> cells;
extern std::unordered_map<int32_t, int32_t> universe_map;
extern vector<unique_ptr<Universe>> universes;
} // namespace model
//==============================================================================
//! A geometry primitive that fills all space and contains cells.
//==============================================================================
class Universe {
public:
int32_t id_; //!< Unique ID
vector<int32_t> cells_; //!< Cells within this universe
//! \brief Write universe information to an HDF5 group.
//! \param group_id An HDF5 group id.
virtual void to_hdf5(hid_t group_id) const;
virtual bool find_cell(Particle& p) const;
BoundingBox bounding_box() const;
const GeometryType& geom_type() const { return geom_type_; }
GeometryType& geom_type() { return geom_type_; }
unique_ptr<UniversePartitioner> partitioner_;
private:
GeometryType geom_type_ = GeometryType::CSG;
};
//==============================================================================
//==============================================================================
@ -296,35 +269,6 @@ protected:
vector<int32_t>::iterator start, const vector<int32_t>& rpn);
};
//==============================================================================
//! Speeds up geometry searches by grouping cells in a search tree.
//
//! Currently this object only works with universes that are divided up by a
//! bunch of z-planes. It could be generalized to other planes, cylinders,
//! and spheres.
//==============================================================================
class UniversePartitioner {
public:
explicit UniversePartitioner(const Universe& univ);
//! Return the list of cells that could contain the given coordinates.
const vector<int32_t>& get_cells(Position r, Direction u) const;
private:
//! A sorted vector of indices to surfaces that partition the universe
vector<int32_t> surfs_;
//! Vectors listing the indices of the cells that lie within each partition
//
//! There are n+1 partitions with n surfaces. `partitions_.front()` gives the
//! cells that lie on the negative side of `surfs_.front()`.
//! `partitions_.back()` gives the cells that lie on the positive side of
//! `surfs_.back()`. Otherwise, `partitions_[i]` gives cells sandwiched
//! between `surfs_[i-1]` and `surfs_[i]`.
vector<vector<int32_t>> partitions_;
};
//==============================================================================
//! Define an instance of a particular cell
//==============================================================================
@ -356,9 +300,5 @@ struct CellInstanceHash {
void read_cells(pugi::xml_node node);
#ifdef DAGMC
class DAGUniverse;
#endif
} // namespace openmc
#endif // OPENMC_CELL_H

78
include/openmc/universe.h Normal file
View file

@ -0,0 +1,78 @@
#ifndef OPENMC_UNIVERSE_H
#define OPENMC_UNIVERSE_H
#include "openmc/cell.h"
namespace openmc {
#ifdef DAGMC
class DAGUniverse;
#endif
class Universe;
class UniversePartitioner;
namespace model {
extern std::unordered_map<int32_t, int32_t> universe_map;
extern vector<unique_ptr<Universe>> universes;
} // namespace model
//==============================================================================
//! A geometry primitive that fills all space and contains cells.
//==============================================================================
class Universe {
public:
int32_t id_; //!< Unique ID
vector<int32_t> cells_; //!< Cells within this universe
//! \brief Write universe information to an HDF5 group.
//! \param group_id An HDF5 group id.
virtual void to_hdf5(hid_t group_id) const;
virtual bool find_cell(Particle& p) const;
BoundingBox bounding_box() const;
const GeometryType& geom_type() const { return geom_type_; }
GeometryType& geom_type() { return geom_type_; }
unique_ptr<UniversePartitioner> partitioner_;
private:
GeometryType geom_type_ = GeometryType::CSG;
};
//==============================================================================
//! Speeds up geometry searches by grouping cells in a search tree.
//
//! Currently this object only works with universes that are divided up by a
//! bunch of z-planes. It could be generalized to other planes, cylinders,
//! and spheres.
//==============================================================================
class UniversePartitioner {
public:
explicit UniversePartitioner(const Universe& univ);
//! Return the list of cells that could contain the given coordinates.
const vector<int32_t>& get_cells(Position r, Direction u) const;
private:
//! A sorted vector of indices to surfaces that partition the universe
vector<int32_t> surfs_;
//! Vectors listing the indices of the cells that lie within each partition
//
//! There are n+1 partitions with n surfaces. `partitions_.front()` gives the
//! cells that lie on the negative side of `surfs_.front()`.
//! `partitions_.back()` gives the cells that lie on the positive side of
//! `surfs_.back()`. Otherwise, `partitions_[i]` gives cells sandwiched
//! between `surfs_[i-1]` and `surfs_[i]`.
vector<vector<int32_t>> partitions_;
};
} // namespace openmc
#endif // OPENMC_UNIVERSE_H

2
openmc/data/decay.py
View file

@ -442,7 +442,7 @@ class Decay(EqualityMixin):
# Read continuous spectrum
ci = {}
params, ci['probability'] = get_tab1_record(file_obj)
ci['type'] = get_decay_modes(params[0])
ci['from_mode'] = get_decay_modes(params[0])
# Read covariance (Ek, Fk) table
LCOV = params[3]

4
openmc/executor.py
View file

@ -166,13 +166,13 @@ def plot_inline(plots, openmc_exec='openmc', cwd='.'):
plots = [plots]
# Create plots.xml
openmc.Plots(plots).export_to_xml()
openmc.Plots(plots).export_to_xml(cwd)
# Run OpenMC in geometry plotting mode
plot_geometry(False, openmc_exec, cwd)
if plots is not None:
images = [_get_plot_image(p) for p in plots]
images = [_get_plot_image(p, cwd) for p in plots]
display(*images)

8
openmc/model/model.py
View file

@ -395,7 +395,7 @@ class Model:
depletion_operator.cleanup_when_done = True
depletion_operator.finalize()
def export_to_xml(self, directory='.'):
def export_to_xml(self, directory='.', remove_surfs=False):
"""Export model to XML files.
Parameters
@ -403,7 +403,11 @@ class Model:
directory : str
Directory to write XML files to. If it doesn't exist already, it
will be created.
remove_surfs : bool
Whether or not to remove redundant surfaces from the geometry when
exporting.
.. versionadded:: 0.13.1
"""
# Create directory if required
d = Path(directory)
@ -411,7 +415,7 @@ class Model:
d.mkdir(parents=True)
self.settings.export_to_xml(d)
self.geometry.export_to_xml(d)
self.geometry.export_to_xml(d, remove_surfs=remove_surfs)
# If a materials collection was specified, export it. Otherwise, look
# for all materials in the geometry and use that to automatically build

139
openmc/model/surface_composite.py
View file

@ -6,7 +6,6 @@ import numpy as np
import openmc
from openmc.checkvalue import check_greater_than, check_value
class CompositeSurface(ABC):
"""Multiple primitive surfaces combined into a composite surface"""
@ -59,11 +58,11 @@ class CylinderSector(CompositeSurface):
This class acts as a proper surface, meaning that unary `+` and `-`
operators applied to it will produce a half-space. The negative
side is defined to be the region inside of the cylinder sector.
Parameters
----------
center : iterable of float
Coordinate for central axes of cylinders in the (y, z), (z, x), or
center : iterable of float
Coordinate for central axes of cylinders in the (y, z), (z, x), or
(x, y) basis. Defaults to (0,0)
r1 : float
Inner cylinder radii
@ -137,6 +136,138 @@ class CylinderSector(CompositeSurface):
def __pos__(self):
return +self.outer_cyl | -self.inner_cyl | +self.plane0 | -self.plane1
class IsogonalOctagon(CompositeSurface):
"""Infinite isogonal octagon composite surface
An isogonal octagon is composed of eight planar surfaces. The prism is
parallel to the x, y, or z axis. The remaining two axes (y and z, z and x,
or x and y) serve as a basis for constructing the surfaces. Two surfaces
are parallel to the first basis axis, two surfaces are parallel
to the second basis axis, and the remaining four surfaces intersect both
basis axes at 45 degree angles.
This class acts as a proper surface, meaning that unary `+` and `-`
operators applied to it will produce a half-space. The negative side is
defined to be the region inside of the octogonal prism.
Parameters
----------
center : iterable of float
Coordinate for the central axis of the octagon in the
(y, z), (z, x), or (x, y) basis.
r1 : float
Half-width of octagon across its basis axis-parallel sides in units
of cm. Must be less than :math:`r_2\sqrt{2}`.
r2 : float
Half-width of octagon across its basis axis intersecting sides in
units of cm. Must be less than than :math:`r_1\sqrt{2}`.
axis : {'x', 'y', 'z'}
Central axis of octagon. Defaults to 'z'
**kwargs
Keyword arguments passed to underlying plane classes
Attributes
----------
top : openmc.ZPlane, openmc.XPlane, or openmc.YPlane
Top planar surface of octagon
bottom : openmc.ZPlane, openmc.XPlane, or openmc.YPlane
Bottom planar surface of octagon
right : openmc.YPlane, openmc.ZPlane, or openmc.XPlane
Right planar surface of octagon
left : openmc.YPlane, openmc.ZPlane, or openmc.XPlane
Left planar surface of octagon
upper_right : openmc.Plane
Upper right planar surface of octagon
lower_right : openmc.Plane
Lower right planar surface of octagon
lower_left : openmc.Plane
Lower left planar surface of octagon
upper_left : openmc.Plane
Upper left planar surface of octagon
"""
_surface_names = ('top', 'bottom',
'upper_right', 'lower_left',
'right', 'left',
'lower_right', 'upper_left')
def __init__(self, center, r1, r2, axis='z', **kwargs):
c1, c2 = center
# Coords for axis-perpendicular planes
ctop = c1 + r1
cbottom = c1 - r1
cright = c2 + r1
cleft = c2 - r1
# Side lengths
if r2 > r1 * sqrt(2):
raise ValueError(f'r2 is greater than sqrt(2) * r1. Octagon' + \
' may be erroneous.')
if r1 > r2 * sqrt(2):
raise ValueError(f'r1 is greater than sqrt(2) * r2. Octagon' + \
' may be erroneous.')
L_basis_ax = (r2 * sqrt(2) - r1)
# Coords for quadrant planes
p1_ur = np.array([L_basis_ax, r1, 0.])
p2_ur = np.array([r1, L_basis_ax, 0.])
p3_ur = np.array([r1, L_basis_ax, 1.])
p1_lr = np.array([r1, -L_basis_ax, 0.])
p2_lr = np.array([L_basis_ax, -r1, 0.])
p3_lr = np.array([L_basis_ax, -r1, 1.])
points = [p1_ur, p2_ur, p3_ur, p1_lr, p2_lr, p3_lr]
# Orientation specific variables
if axis == 'z':
coord_map = [0, 1, 2]
self.top = openmc.YPlane(ctop, **kwargs)
self.bottom = openmc.YPlane(cbottom, **kwargs)
self.right = openmc.XPlane(cright, **kwargs)
self.left = openmc.XPlane(cleft, **kwargs)
elif axis == 'y':
coord_map = [1, 2, 0]
self.top = openmc.XPlane(ctop, **kwargs)
self.bottom = openmc.XPlane(cbottom, **kwargs)
self.right = openmc.ZPlane(cright, **kwargs)
self.left = openmc.ZPlane(cleft, **kwargs)
elif axis == 'x':
coord_map = [2, 0, 1]
self.top = openmc.ZPlane(ctop, **kwargs)
self.bottom = openmc.ZPlane(cbottom, **kwargs)
self.right = openmc.YPlane(cright, **kwargs)
self.left = openmc.YPlane(cleft, **kwargs)
# Put our coordinates in (x,y,z) order
for p in points:
p[:] = p[coord_map]
self.upper_right = openmc.Plane.from_points(p1_ur, p2_ur, p3_ur,
**kwargs)
self.lower_right = openmc.Plane.from_points(p1_lr, p2_lr, p3_lr,
**kwargs)
self.lower_left = openmc.Plane.from_points(-p1_ur, -p2_ur, -p3_ur,
**kwargs)
self.upper_left = openmc.Plane.from_points(-p1_lr, -p2_lr, -p3_lr,
**kwargs)
def __neg__(self):
return -self.top & +self.bottom & -self.right & +self.left & \
+self.upper_right & +self.lower_right & -self.lower_left & \
-self.upper_left
def __pos__(self):
return +self.top | -self.bottom | +self.right | -self.left | \
-self.upper_right | -self.lower_right | +self.lower_left | \
+self.upper_left
class RightCircularCylinder(CompositeSurface):

12
openmc/plots.py
View file

@ -164,16 +164,16 @@ _SVG_COLORS = {
}
def _get_plot_image(plot):
def _get_plot_image(plot, cwd):
from IPython.display import Image
# Make sure .png file was created
stem = plot.filename if plot.filename is not None else f'plot_{plot.id}'
png_file = f'{stem}.png'
if not Path(png_file).exists():
png_file = Path(cwd) / f'{stem}.png'
if not png_file.exists():
raise FileNotFoundError(f"Could not find .png image for plot {plot.id}")
return Image(png_file)
return Image(str(png_file))
class Plot(IDManagerMixin):
@ -784,13 +784,13 @@ class Plot(IDManagerMixin):
"""
# Create plots.xml
Plots([self]).export_to_xml()
Plots([self]).export_to_xml(cwd)
# Run OpenMC in geometry plotting mode
openmc.plot_geometry(False, openmc_exec, cwd)
# Return produced image
return _get_plot_image(self)
return _get_plot_image(self, cwd)
class Plots(cv.CheckedList):

130
openmc/universe.py
View file

@ -1,9 +1,10 @@
from abc import ABC, abstractmethod
from collections import OrderedDict
from collections.abc import Iterable
from copy import copy, deepcopy
from copy import deepcopy
from numbers import Real
import random
from pathlib import Path
from tempfile import TemporaryDirectory
from xml.etree import ElementTree as ET
import numpy as np
@ -267,13 +268,9 @@ class Universe(UniverseBase):
def plot(self, origin=(0., 0., 0.), width=(1., 1.), pixels=(200, 200),
basis='xy', color_by='cell', colors=None, seed=None,
**kwargs):
openmc_exec='openmc', **kwargs):
"""Display a slice plot of the universe.
To display or save the plot, call :func:`matplotlib.pyplot.show` or
:func:`matplotlib.pyplot.savefig`. In a Jupyter notebook, enabling the
matplotlib inline backend will show the plot inline.
Parameters
----------
origin : Iterable of float
@ -297,14 +294,14 @@ class Universe(UniverseBase):
# Make water blue
water = openmc.Cell(fill=h2o)
universe.plot(..., colors={water: (0., 0., 1.))
seed : int
Seed for the random number generator
openmc_exec : str
Path to OpenMC executable.
seed : hashable object or None
Hashable object which is used to seed the random number generator
used to select colors. If None, the generator is seeded from the
current time.
.. versionadded:: 0.13.1
**kwargs
All keyword arguments are passed to
:func:`matplotlib.pyplot.imshow`.
Keyword arguments passed to :func:`matplotlib.pyplot.imshow`
Returns
-------
@ -312,83 +309,54 @@ class Universe(UniverseBase):
Resulting image
"""
import matplotlib.image as mpimg
import matplotlib.pyplot as plt
# Seed the random number generator
if seed is not None:
random.seed(seed)
if colors is None:
# Create default dictionary if none supplied
colors = {}
else:
# Convert to RGBA if necessary
colors = copy(colors)
for obj, color in colors.items():
if isinstance(color, str):
if color.lower() not in _SVG_COLORS:
raise ValueError(f"'{color}' is not a valid color.")
colors[obj] = [x/255 for x in
_SVG_COLORS[color.lower()]] + [1.0]
elif len(color) == 3:
colors[obj] = list(color) + [1.0]
# Determine extents of plot
if basis == 'xy':
x_min = origin[0] - 0.5*width[0]
x_max = origin[0] + 0.5*width[0]
y_min = origin[1] - 0.5*width[1]
y_max = origin[1] + 0.5*width[1]
x, y = 0, 1
elif basis == 'yz':
# The x-axis will correspond to physical y and the y-axis will
# correspond to physical z
x_min = origin[1] - 0.5*width[0]
x_max = origin[1] + 0.5*width[0]
y_min = origin[2] - 0.5*width[1]
y_max = origin[2] + 0.5*width[1]
x, y = 1, 2
elif basis == 'xz':
# The y-axis will correspond to physical z
x_min = origin[0] - 0.5*width[0]
x_max = origin[0] + 0.5*width[0]
y_min = origin[2] - 0.5*width[1]
y_max = origin[2] + 0.5*width[1]
x, y = 0, 2
x_min = origin[x] - 0.5*width[0]
x_max = origin[x] + 0.5*width[0]
y_min = origin[y] - 0.5*width[1]
y_max = origin[y] + 0.5*width[1]
# Determine locations to determine cells at
x_coords = np.linspace(x_min, x_max, pixels[0], endpoint=False) + \
0.5*(x_max - x_min)/pixels[0]
y_coords = np.linspace(y_max, y_min, pixels[1], endpoint=False) - \
0.5*(y_max - y_min)/pixels[1]
with TemporaryDirectory() as tmpdir:
model = openmc.Model()
model.geometry = openmc.Geometry(self)
if seed is not None:
model.settings.seed = seed
# Initialize output image in RGBA format. Flip the pixels from
# traditional (x, y) to (y, x) used in graphics.
img = np.zeros((pixels[1], pixels[0], 4))
for i, x in enumerate(x_coords):
for j, y in enumerate(y_coords):
if basis == 'xy':
path = self.find((x, y, origin[2]))
elif basis == 'yz':
path = self.find((origin[0], x, y))
elif basis == 'xz':
path = self.find((x, origin[1], y))
# Create plot object matching passed arguments
plot = openmc.Plot()
plot.origin = origin
plot.width = width
plot.pixels = pixels
plot.basis = basis
plot.color_by = color_by
plot.colors = colors
model.plots.append(plot)
if len(path) > 0:
try:
if color_by == 'cell':
obj = path[-1]
elif color_by == 'material':
if path[-1].fill_type == 'material':
obj = path[-1].fill
else:
continue
except AttributeError:
continue
if obj not in colors:
colors[obj] = (random.random(), random.random(),
random.random(), 1.0)
img[j, i, :] = colors[obj]
# Run OpenMC in geometry plotting mode
model.plot_geometry(False, cwd=tmpdir, openmc_exec=openmc_exec)
# Display image
return plt.imshow(img, extent=(x_min, x_max, y_min, y_max),
interpolation='nearest', **kwargs)
# Read image from file
img = mpimg.imread(Path(tmpdir) / f'plot_{plot.id}.png')
# Create a figure sized such that the size of the axes within
# exactly matches the number of pixels specified
px = 1/plt.rcParams['figure.dpi']
fig, ax = plt.subplots()
params = fig.subplotpars
width = pixels[0]*px/(params.right - params.left)
height = pixels[0]*px/(params.top - params.bottom)
fig.set_size_inches(width, height)
# Plot image and return the axes
return ax.imshow(img, extent=(x_min, x_max, y_min, y_max), **kwargs)
def add_cell(self, cell):
"""Add a cell to the universe.

207
src/cell.cpp
View file

@ -34,8 +34,6 @@ namespace model {
std::unordered_map<int32_t, int32_t> cell_map;
vector<unique_ptr<Cell>> cells;
std::unordered_map<int32_t, int32_t> universe_map;
vector<unique_ptr<Universe>> universes;
} // namespace model
//==============================================================================
@ -186,66 +184,6 @@ vector<int32_t> generate_rpn(int32_t cell_id, vector<int32_t> infix)
return rpn;
}
//==============================================================================
// Universe implementation
//==============================================================================
void Universe::to_hdf5(hid_t universes_group) const
{
// Create a group for this universe.
auto group = create_group(universes_group, fmt::format("universe {}", id_));
// Write the geometry representation type.
write_string(group, "geom_type", "csg", false);
// Write the contained cells.
if (cells_.size() > 0) {
vector<int32_t> cell_ids;
for (auto i_cell : cells_)
cell_ids.push_back(model::cells[i_cell]->id_);
write_dataset(group, "cells", cell_ids);
}
close_group(group);
}
bool Universe::find_cell(Particle& p) const
{
const auto& cells {
!partitioner_ ? cells_ : partitioner_->get_cells(p.r_local(), p.u_local())};
for (auto it = cells.begin(); it != cells.end(); it++) {
int32_t i_cell = *it;
int32_t i_univ = p.coord(p.n_coord() - 1).universe;
if (model::cells[i_cell]->universe_ != i_univ)
continue;
// Check if this cell contains the particle;
Position r {p.r_local()};
Direction u {p.u_local()};
auto surf = p.surface();
if (model::cells[i_cell]->contains(r, u, surf)) {
p.coord(p.n_coord() - 1).cell = i_cell;
return true;
}
}
return false;
}
BoundingBox Universe::bounding_box() const
{
BoundingBox bbox = {INFTY, -INFTY, INFTY, -INFTY, INFTY, -INFTY};
if (cells_.size() == 0) {
return {};
} else {
for (const auto& cell : cells_) {
auto& c = model::cells[cell];
bbox |= c->bounding_box();
}
}
return bbox;
}
//==============================================================================
// Cell implementation
//==============================================================================
@ -876,151 +814,6 @@ bool CSGCell::contains_complex(
}
}
//==============================================================================
// UniversePartitioner implementation
//==============================================================================
UniversePartitioner::UniversePartitioner(const Universe& univ)
{
// Define an ordered set of surface indices that point to z-planes. Use a
// functor to to order the set by the z0_ values of the corresponding planes.
struct compare_surfs {
bool operator()(const int32_t& i_surf, const int32_t& j_surf) const
{
const auto* surf = model::surfaces[i_surf].get();
const auto* zplane = dynamic_cast<const SurfaceZPlane*>(surf);
double zi = zplane->z0_;
surf = model::surfaces[j_surf].get();
zplane = dynamic_cast<const SurfaceZPlane*>(surf);
double zj = zplane->z0_;
return zi < zj;
}
};
std::set<int32_t, compare_surfs> surf_set;
// Find all of the z-planes in this universe. A set is used here for the
// O(log(n)) insertions that will ensure entries are not repeated.
for (auto i_cell : univ.cells_) {
for (auto token : model::cells[i_cell]->rpn_) {
if (token < OP_UNION) {
auto i_surf = std::abs(token) - 1;
const auto* surf = model::surfaces[i_surf].get();
if (const auto* zplane = dynamic_cast<const SurfaceZPlane*>(surf))
surf_set.insert(i_surf);
}
}
}
// Populate the surfs_ vector from the ordered set.
surfs_.insert(surfs_.begin(), surf_set.begin(), surf_set.end());
// Populate the partition lists.
partitions_.resize(surfs_.size() + 1);
for (auto i_cell : univ.cells_) {
// It is difficult to determine the bounds of a complex cell, so add complex
// cells to all partitions.
if (!model::cells[i_cell]->simple_) {
for (auto& p : partitions_)
p.push_back(i_cell);
continue;
}
// Find the tokens for bounding z-planes.
int32_t lower_token = 0, upper_token = 0;
double min_z, max_z;
for (auto token : model::cells[i_cell]->rpn_) {
if (token < OP_UNION) {
const auto* surf = model::surfaces[std::abs(token) - 1].get();
if (const auto* zplane = dynamic_cast<const SurfaceZPlane*>(surf)) {
if (lower_token == 0 || zplane->z0_ < min_z) {
lower_token = token;
min_z = zplane->z0_;
}
if (upper_token == 0 || zplane->z0_ > max_z) {
upper_token = token;
max_z = zplane->z0_;
}
}
}
}
// If there are no bounding z-planes, add this cell to all partitions.
if (lower_token == 0) {
for (auto& p : partitions_)
p.push_back(i_cell);
continue;
}
// Find the first partition this cell lies in. If the lower_token indicates
// a negative halfspace, then the cell is unbounded in the lower direction
// and it lies in the first partition onward. Otherwise, it is bounded by
// the positive halfspace given by the lower_token.
int first_partition = 0;
if (lower_token > 0) {
for (int i = 0; i < surfs_.size(); ++i) {
if (lower_token == surfs_[i] + 1) {
first_partition = i + 1;
break;
}
}
}
// Find the last partition this cell lies in. The logic is analogous to the
// logic for first_partition.
int last_partition = surfs_.size();
if (upper_token < 0) {
for (int i = first_partition; i < surfs_.size(); ++i) {
if (upper_token == -(surfs_[i] + 1)) {
last_partition = i;
break;
}
}
}
// Add the cell to all relevant partitions.
for (int i = first_partition; i <= last_partition; ++i) {
partitions_[i].push_back(i_cell);
}
}
}
const vector<int32_t>& UniversePartitioner::get_cells(
Position r, Direction u) const
{
// Perform a binary search for the partition containing the given coordinates.
int left = 0;
int middle = (surfs_.size() - 1) / 2;
int right = surfs_.size() - 1;
while (true) {
// Check the sense of the coordinates for the current surface.
const auto& surf = *model::surfaces[surfs_[middle]];
if (surf.sense(r, u)) {
// The coordinates lie in the positive halfspace. Recurse if there are
// more surfaces to check. Otherwise, return the cells on the positive
// side of this surface.
int right_leaf = right - (right - middle) / 2;
if (right_leaf != middle) {
left = middle + 1;
middle = right_leaf;
} else {
return partitions_[middle + 1];
}
} else {
// The coordinates lie in the negative halfspace. Recurse if there are
// more surfaces to check. Otherwise, return the cells on the negative
// side of this surface.
int left_leaf = left + (middle - left) / 2;
if (left_leaf != middle) {
right = middle - 1;
middle = left_leaf;
} else {
return partitions_[middle];
}
}
}
}
//==============================================================================
// Non-method functions
//==============================================================================

221
src/universe.cpp Normal file
View file

@ -0,0 +1,221 @@
#include "openmc/universe.h"
#include <set>
#include "openmc/hdf5_interface.h"
namespace openmc {
namespace model {
std::unordered_map<int32_t, int32_t> universe_map;
vector<unique_ptr<Universe>> universes;
} // namespace model
//==============================================================================
// Universe implementation
//==============================================================================
void Universe::to_hdf5(hid_t universes_group) const
{
// Create a group for this universe.
auto group = create_group(universes_group, fmt::format("universe {}", id_));
// Write the geometry representation type.
write_string(group, "geom_type", "csg", false);
// Write the contained cells.
if (cells_.size() > 0) {
vector<int32_t> cell_ids;
for (auto i_cell : cells_)
cell_ids.push_back(model::cells[i_cell]->id_);
write_dataset(group, "cells", cell_ids);
}
close_group(group);
}
bool Universe::find_cell(Particle& p) const
{
const auto& cells {
!partitioner_ ? cells_ : partitioner_->get_cells(p.r_local(), p.u_local())};
for (auto it = cells.begin(); it != cells.end(); it++) {
int32_t i_cell = *it;
int32_t i_univ = p.coord(p.n_coord() - 1).universe;
if (model::cells[i_cell]->universe_ != i_univ)
continue;
// Check if this cell contains the particle;
Position r {p.r_local()};
Direction u {p.u_local()};
auto surf = p.surface();
if (model::cells[i_cell]->contains(r, u, surf)) {
p.coord(p.n_coord() - 1).cell = i_cell;
return true;
}
}
return false;
}
BoundingBox Universe::bounding_box() const
{
BoundingBox bbox = {INFTY, -INFTY, INFTY, -INFTY, INFTY, -INFTY};
if (cells_.size() == 0) {
return {};
} else {
for (const auto& cell : cells_) {
auto& c = model::cells[cell];
bbox |= c->bounding_box();
}
}
return bbox;
}
//==============================================================================
// UniversePartitioner implementation
//==============================================================================
UniversePartitioner::UniversePartitioner(const Universe& univ)
{
// Define an ordered set of surface indices that point to z-planes. Use a
// functor to to order the set by the z0_ values of the corresponding planes.
struct compare_surfs {
bool operator()(const int32_t& i_surf, const int32_t& j_surf) const
{
const auto* surf = model::surfaces[i_surf].get();
const auto* zplane = dynamic_cast<const SurfaceZPlane*>(surf);
double zi = zplane->z0_;
surf = model::surfaces[j_surf].get();
zplane = dynamic_cast<const SurfaceZPlane*>(surf);
double zj = zplane->z0_;
return zi < zj;
}
};
std::set<int32_t, compare_surfs> surf_set;
// Find all of the z-planes in this universe. A set is used here for the
// O(log(n)) insertions that will ensure entries are not repeated.
for (auto i_cell : univ.cells_) {
for (auto token : model::cells[i_cell]->rpn_) {
if (token < OP_UNION) {
auto i_surf = std::abs(token) - 1;
const auto* surf = model::surfaces[i_surf].get();
if (const auto* zplane = dynamic_cast<const SurfaceZPlane*>(surf))
surf_set.insert(i_surf);
}
}
}
// Populate the surfs_ vector from the ordered set.
surfs_.insert(surfs_.begin(), surf_set.begin(), surf_set.end());
// Populate the partition lists.
partitions_.resize(surfs_.size() + 1);
for (auto i_cell : univ.cells_) {
// It is difficult to determine the bounds of a complex cell, so add complex
// cells to all partitions.
if (!model::cells[i_cell]->simple_) {
for (auto& p : partitions_)
p.push_back(i_cell);
continue;
}
// Find the tokens for bounding z-planes.
int32_t lower_token = 0, upper_token = 0;
double min_z, max_z;
for (auto token : model::cells[i_cell]->rpn_) {
if (token < OP_UNION) {
const auto* surf = model::surfaces[std::abs(token) - 1].get();
if (const auto* zplane = dynamic_cast<const SurfaceZPlane*>(surf)) {
if (lower_token == 0 || zplane->z0_ < min_z) {
lower_token = token;
min_z = zplane->z0_;
}
if (upper_token == 0 || zplane->z0_ > max_z) {
upper_token = token;
max_z = zplane->z0_;
}
}
}
}
// If there are no bounding z-planes, add this cell to all partitions.
if (lower_token == 0) {
for (auto& p : partitions_)
p.push_back(i_cell);
continue;
}
// Find the first partition this cell lies in. If the lower_token indicates
// a negative halfspace, then the cell is unbounded in the lower direction
// and it lies in the first partition onward. Otherwise, it is bounded by
// the positive halfspace given by the lower_token.
int first_partition = 0;
if (lower_token > 0) {
for (int i = 0; i < surfs_.size(); ++i) {
if (lower_token == surfs_[i] + 1) {
first_partition = i + 1;
break;
}
}
}
// Find the last partition this cell lies in. The logic is analogous to the
// logic for first_partition.
int last_partition = surfs_.size();
if (upper_token < 0) {
for (int i = first_partition; i < surfs_.size(); ++i) {
if (upper_token == -(surfs_[i] + 1)) {
last_partition = i;
break;
}
}
}
// Add the cell to all relevant partitions.
for (int i = first_partition; i <= last_partition; ++i) {
partitions_[i].push_back(i_cell);
}
}
}
const vector<int32_t>& UniversePartitioner::get_cells(
Position r, Direction u) const
{
// Perform a binary search for the partition containing the given coordinates.
int left = 0;
int middle = (surfs_.size() - 1) / 2;
int right = surfs_.size() - 1;
while (true) {
// Check the sense of the coordinates for the current surface.
const auto& surf = *model::surfaces[surfs_[middle]];
if (surf.sense(r, u)) {
// The coordinates lie in the positive halfspace. Recurse if there are
// more surfaces to check. Otherwise, return the cells on the positive
// side of this surface.
int right_leaf = right - (right - middle) / 2;
if (right_leaf != middle) {
left = middle + 1;
middle = right_leaf;
} else {
return partitions_[middle + 1];
}
} else {
// The coordinates lie in the negative halfspace. Recurse if there are
// more surfaces to check. Otherwise, return the cells on the negative
// side of this surface.
int left_leaf = left + (middle - left) / 2;
if (left_leaf != middle) {
right = middle - 1;
middle = left_leaf;
} else {
return partitions_[middle];
}
}
}
}
} // namespace openmc

6
tests/regression_tests/cpp_driver/test.py
View file

@ -33,12 +33,14 @@ def cpp_driver(request):
os.chdir(str(local_builddir))
if config['mpi']:
os.environ['CXX'] = 'mpicxx'
mpi_arg = "On"
else:
mpi_arg = "Off"
try:
print("Building driver")
# Run cmake/make to build the shared libary
subprocess.run(['cmake', os.path.pardir], check=True)
subprocess.run(['cmake', os.path.pardir, f'-DOPENMC_USE_MPI={mpi_arg}'], check=True)
subprocess.run(['make'], check=True)
os.chdir(os.path.pardir)

6
tests/regression_tests/external_moab/test.py
View file

@ -49,12 +49,14 @@ def cpp_driver(request):
os.chdir(str(local_builddir))
if config['mpi']:
os.environ['CXX'] = 'mpicxx'
mpi_arg = "On"
else:
mpi_arg = "Off"
try:
print("Building driver")
# Run cmake/make to build the shared libary
subprocess.run(['cmake', os.path.pardir], check=True)
subprocess.run(['cmake', os.path.pardir, f'-DOPENMC_USE_MPI={mpi_arg}'], check=True)
subprocess.run(['make'], check=True)
os.chdir(os.path.pardir)

6
tests/unit_tests/test_lattice.py
View file

@ -157,11 +157,11 @@ def hlat3(pincell1, pincell2, uo2, water, zr):
def test_get_nuclides(rlat2, rlat3, hlat2, hlat3):
for lat in (rlat2, hlat2):
nucs = rlat2.get_nuclides()
nucs = lat.get_nuclides()
assert sorted(nucs) == ['H1', 'O16', 'U235',
'Zr90', 'Zr91', 'Zr92', 'Zr94', 'Zr96']
for lat in (rlat3, hlat3):
nucs = rlat3.get_nuclides()
nucs = lat.get_nuclides()
assert sorted(nucs) == ['H1', 'H2', 'O16', 'U235',
'Zr90', 'Zr91', 'Zr92', 'Zr94', 'Zr96']
@ -360,7 +360,7 @@ def test_show_indices():
lines = openmc.HexLattice.show_indices(i).split('\n')
assert len(lines) == 4*i - 3
lines_x = openmc.HexLattice.show_indices(i, 'x').split('\n')
assert len(lines) == 4*i - 3
assert len(lines_x) == 4*i - 3
def test_unset_universes():

84
tests/unit_tests/test_surface_composite.py
View file

@ -151,6 +151,7 @@ def test_cone_one_sided(axis, point_pos, point_neg, ll_true):
# Make sure repr works
repr(s)
@pytest.mark.parametrize(
"axis, indices", [
("X", [0, 1, 2]),
@ -170,7 +171,7 @@ def test_cylinder_sector(axis, indices):
assert isinstance(s.inner_cyl, getattr(openmc, axis + "Cylinder"))
assert isinstance(s.plane0, openmc.Plane)
assert isinstance(s.plane1, openmc.Plane)
# Make sure boundary condition propagates
s.boundary_type = 'reflective'
assert s.boundary_type == 'reflective'
@ -178,7 +179,7 @@ def test_cylinder_sector(axis, indices):
assert s.inner_cyl.boundary_type == 'reflective'
assert s.plane0.boundary_type == 'reflective'
assert s.plane1.boundary_type == 'reflective'
# Check bounding box
ll, ur = (+s).bounding_box
assert np.all(np.isinf(ll))
@ -202,7 +203,86 @@ def test_cylinder_sector(axis, indices):
assert ur_t == pytest.approx(ur + t)
assert ll_t == pytest.approx(ll + t)
# Check invalid r1, r2 combinations
with pytest.raises(ValueError):
openmc.model.IsogonalOctagon(center, r1=1.0, r2=10.)
with pytest.raises(ValueError):
openmc.model.IsogonalOctagon(center, r1=10., r2=1.)
# Make sure repr works
repr(s)
@pytest.mark.parametrize(
"axis, plane_tb, plane_lr, axis_idx", [
("x", "Z", "Y", 0),
("y", "X", "Z", 1),
("z", "Y", "X", 2),
]
)
def test_isogonal_octagon(axis, plane_tb, plane_lr, axis_idx):
center = np.array([0., 0.])
point_pos = np.array([0.8, 0.8])
point_neg = np.array([0.7, 0.7])
r1 = 1.
r2 = 1.
plane_top_bottom = getattr(openmc, plane_tb + "Plane")
plane_left_right = getattr(openmc, plane_lr + "Plane")
s = openmc.model.IsogonalOctagon(center, r1, r2, axis=axis)
assert isinstance(s.top, plane_top_bottom)
assert isinstance(s.bottom, plane_top_bottom)
assert isinstance(s.right, plane_left_right)
assert isinstance(s.left, plane_left_right)
assert isinstance(s.upper_right, openmc.Plane)
assert isinstance(s.lower_right, openmc.Plane)
assert isinstance(s.upper_left, openmc.Plane)
assert isinstance(s.lower_left, openmc.Plane)
# Make sure boundary condition propagates
s.boundary_type = 'reflective'
assert s.boundary_type == 'reflective'
assert s.top.boundary_type == 'reflective'
assert s.bottom.boundary_type == 'reflective'
assert s.right.boundary_type == 'reflective'
assert s.left.boundary_type == 'reflective'
assert s.upper_right.boundary_type == 'reflective'
assert s.lower_right.boundary_type == 'reflective'
assert s.lower_left.boundary_type == 'reflective'
assert s.upper_left.boundary_type == 'reflective'
# Check bounding box
center = np.insert(center, axis_idx, np.inf)
xmax, ymax, zmax = center + r1
coord_min = center - r1
coord_min[axis_idx] *= -1
xmin, ymin, zmin = coord_min
ll, ur = (+s).bounding_box
assert np.all(np.isinf(ll))
assert np.all(np.isinf(ur))
ll, ur = (-s).bounding_box
assert ur == pytest.approx((xmax, ymax, zmax))
assert ll == pytest.approx((xmin, ymin, zmin))
# __contains__ on associated half-spaces
point_pos = np.insert(point_pos, axis_idx, 0)
point_neg = np.insert(point_neg, axis_idx, 0)
assert point_pos in +s
assert point_pos not in -s
assert point_neg in -s
assert point_neg not in +s
# translate method
t = uniform(-5.0, 5.0)
s_t = s.translate((t, t, t))
ll_t, ur_t = (-s_t).bounding_box
assert ur_t == pytest.approx(ur + t)
assert ll_t == pytest.approx(ll + t)
# Check invalid r1, r2 combinations
with pytest.raises(ValueError):
openmc.model.IsogonalOctagon(center, r1=1.0, r2=10.)
with pytest.raises(ValueError):
openmc.model.IsogonalOctagon(center, r1=10., r2=1.)
# Make sure repr works
repr(s)

2
tools/ci/gha-install.py
View file

@ -34,7 +34,7 @@ def install(omp=False, mpi=False, phdf5=False, dagmc=False, libmesh=False):
# Use MPI wrappers when building in parallel
if mpi:
os.environ['CXX'] = 'mpicxx'
cmake_cmd.append('-DOPENMC_USE_MPI=on')
# Tell CMake to prefer parallel HDF5 if specified
if phdf5: