diff --git a/CMakeLists.txt b/CMakeLists.txt index ec6b500b4..29eed4fb0 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -9,6 +9,11 @@ set(CMAKE_RUNTIME_OUTPUT_DIRECTORY ${CMAKE_BINARY_DIR}/bin) # Set module path set(CMAKE_MODULE_PATH ${CMAKE_CURRENT_SOURCE_DIR}/cmake/Modules) +# Allow user to specify _ROOT variables +if (NOT (CMAKE_VERSION VERSION_LESS 3.12)) + cmake_policy(SET CMP0074 NEW) +endif() + #=============================================================================== # Command line options #=============================================================================== @@ -35,18 +40,18 @@ endif() #=============================================================================== if(dagmc) find_package(DAGMC REQUIRED) + add_library(dagmc-imported INTERFACE IMPORTED) link_directories(${DAGMC_LIBRARY_DIRS}) + target_link_libraries(dagmc-imported INTERFACE ${DAGMC_LIBRARIES}) + target_include_directories(dagmc-imported INTERFACE ${DAGMC_INCLUDE_DIRS}) +else() + set(DAGMC_FOUND false) endif() #=============================================================================== # HDF5 for binary output #=============================================================================== -# Allow user to specify HDF5_ROOT -if (NOT (CMAKE_VERSION VERSION_LESS 3.12)) - cmake_policy(SET CMP0074 NEW) -endif() - # Unfortunately FindHDF5.cmake will always prefer a serial HDF5 installation # over a parallel installation if both appear on the user's PATH. To get around # this, we check for the environment variable HDF5_ROOT and if it exists, use it @@ -357,8 +362,7 @@ target_link_libraries(libopenmc ${ldflags} ${HDF5_LIBRARIES} ${HDF5_HL_LIBRARIES if(dagmc) target_compile_definitions(libopenmc PRIVATE DAGMC) - target_link_libraries(libopenmc ${DAGMC_LIBRARIES}) - target_include_directories(libopenmc PRIVATE ${DAGMC_INCLUDE_DIRS}) + target_link_libraries(libopenmc dagmc-imported) endif() #=============================================================================== @@ -388,6 +392,8 @@ add_custom_command(TARGET libopenmc POST_BUILD # Install executable, scripts, manpage, license #=============================================================================== +configure_file(cmake/OpenMCConfig.cmake.in "${CMAKE_BINARY_DIR}${CMAKE_FILES_DIRECTORY}/OpenMCConfig.cmake" @ONLY) + set(INSTALL_CONFIGDIR ${CMAKE_INSTALL_LIBDIR}/cmake/OpenMC) install(TARGETS openmc libopenmc faddeeva EXPORT openmc-targets @@ -401,7 +407,7 @@ install(EXPORT openmc-targets DESTINATION ${INSTALL_CONFIGDIR}) install(DIRECTORY src/relaxng DESTINATION ${CMAKE_INSTALL_DATADIR}/openmc) -install(FILES cmake/OpenMCConfig.cmake DESTINATION ${INSTALL_CONFIGDIR}) +install(FILES "${CMAKE_BINARY_DIR}${CMAKE_FILES_DIRECTORY}/OpenMCConfig.cmake" DESTINATION ${INSTALL_CONFIGDIR}) install(FILES man/man1/openmc.1 DESTINATION ${CMAKE_INSTALL_MANDIR}/man1) install(FILES LICENSE DESTINATION "${CMAKE_INSTALL_DOCDIR}" RENAME copyright) install(DIRECTORY include/ DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}) diff --git a/cmake/Modules/FindDAGMC.cmake b/cmake/Modules/FindDAGMC.cmake index c644c8886..bb9bc56dd 100644 --- a/cmake/Modules/FindDAGMC.cmake +++ b/cmake/Modules/FindDAGMC.cmake @@ -12,7 +12,9 @@ find_path(DAGMC_CMAKE_CONFIG NAMES DAGMCConfig.cmake PATHS ENV LD_LIBRARY_PATH PATH_SUFFIXES lib Lib cmake lib/cmake NO_DEFAULT_PATH) - -message(STATUS "Found DAGMC in ${DAGMC_CMAKE_CONFIG}") - -include(${DAGMC_CMAKE_CONFIG}/DAGMCConfig.cmake) +if(DAGMC_CMAKE_CONFIG) + message(STATUS "Found DAGMC in ${DAGMC_CMAKE_CONFIG}") + include(${DAGMC_CMAKE_CONFIG}/DAGMCConfig.cmake) +else() + message(WARNING "Cound not find DAGMC") +endif() diff --git a/cmake/OpenMCConfig.cmake b/cmake/OpenMCConfig.cmake.in similarity index 56% rename from cmake/OpenMCConfig.cmake rename to cmake/OpenMCConfig.cmake.in index 29a0e4542..1c6a488f5 100644 --- a/cmake/OpenMCConfig.cmake +++ b/cmake/OpenMCConfig.cmake.in @@ -5,6 +5,15 @@ find_package(gsl-lite REQUIRED HINTS ${OpenMC_CMAKE_DIR}/../gsl-lite) find_package(pugixml REQUIRED HINTS ${OpenMC_CMAKE_DIR}/../pugixml) find_package(xtl REQUIRED HINTS ${OpenMC_CMAKE_DIR}/../xtl) find_package(xtensor REQUIRED HINTS ${OpenMC_CMAKE_DIR}/../xtensor) +if(@DAGMC_FOUND@) + find_package(DAGMC REQUIRED HINTS @DAGMC_LIBRARY_DIRS@) + add_library(dagmc-imported INTERFACE IMPORTED) + target_link_libraries(dagmc-imported INTERFACE ${DAGMC_LIBRARIES}) + foreach(dir ${DAGMC_LIBRARY_DIRS}) + target_link_libraries(dagmc-imported INTERFACE "-L${dir}") + endforeach() + target_include_directories(dagmc-imported INTERFACE ${DAGMC_INCLUDE_DIRS}) +endif() if(NOT TARGET OpenMC::libopenmc) include("${OpenMC_CMAKE_DIR}/OpenMCTargets.cmake") diff --git a/docs/source/io_formats/settings.rst b/docs/source/io_formats/settings.rst index 2e2648237..83aa0691f 100644 --- a/docs/source/io_formats/settings.rst +++ b/docs/source/io_formats/settings.rst @@ -142,7 +142,7 @@ materials in the problem and is specified using a :ref:`mesh_element`. ---------------------------- Determines whether to use event-based parallelism instead of the default -history-based parallelism. +history-based parallelism. *Default*: false @@ -459,6 +459,15 @@ attributes/sub-elements: *Default*: None + :library: + If this attribute is given, it indicates that the source is to be + instantiated from an externally compiled source function. This source can be + as complex as is required to define the source for your problem. The only + requirement is that there is a function called ``sample_source()``. More + documentation on how to build sources can be found in :ref:`custom_source`. + + *Default*: None + :space: An element specifying the spatial distribution of source sites. This element has the following attributes: @@ -591,6 +600,67 @@ attributes/sub-elements: *Default*: false +.. _custom_source: + +Custom Sources +++++++++++++++ + +It is often the case that one may wish to simulate a complex source +distribution, which may include physics not present within OpenMC or to be phase +space complex. It is possible to define a complex source with an externally +defined source function that is loaded at runtime. A simple example source is +shown below. + +.. code-block:: c++ + + #include "openmc/random_lcg.h" + #include "openmc/source.h" + #include "openmc/particle.h" + + // you must have external C linkage here + extern "C" openmc::Particle::Bank sample_source(uint64_t* seed) { + openmc::Particle::Bank particle; + // weight + particle.particle = openmc::Particle::Type::neutron; + particle.wgt = 1.0; + // position + double angle = 2.0 * M_PI * openmc::prn(seed); + double radius = 3.0; + particle.r.x = radius * std::cos(angle); + particle.r.y = radius * std::sin(angle); + particle.r.z = 0.0; + // angle + particle.u = {1.0, 0.0, 0.0}; + particle.E = 14.08e6; + particle.delayed_group = 0; + return particle; + } + +The above source, creates 14.08 MeV neutrons, with an istropic direction +vector but distributed in a ring with a 3 cm radius. This routine is +not particularly complex, but should serve as an example upon which to build +more complicated sources. + + .. note:: The function signature must be declared to be extern "C". + + .. note:: You should only use the openmc::prn() random number generator + +In order to build your external source, you will need to link it against the +OpenMC shared library. This can be done by writing a CMakeLists.txt file: + +.. code-block:: cmake + + cmake_minimum_required(VERSION 3.3 FATAL_ERROR) + project(openmc_sources CXX) + add_library(source SHARED source_ring.cpp) + find_package(OpenMC REQUIRED HINTS ) + target_link_libraries(source OpenMC::libopenmc) + +After running ``cmake`` and ``make``, you will have a libsource.so (or .dylib) +file in your build directory. Setting the :attr:`openmc.Source.library` +attribute to the path of this shared library will indicate that it should be +used for sampling source particles at runtime. + .. _univariate: Univariate Probability Distributions diff --git a/docs/source/usersguide/install.rst b/docs/source/usersguide/install.rst index 66440a3c3..ec10fea4f 100644 --- a/docs/source/usersguide/install.rst +++ b/docs/source/usersguide/install.rst @@ -404,7 +404,7 @@ to install the Python package in :ref:`"editable" mode `. Prerequisites ------------- -The Python API works with Python 3.4+. In addition to Python itself, the API +The Python API works with Python 3.5+. In addition to Python itself, the API relies on a number of third-party packages. All prerequisites can be installed using Conda_ (recommended), pip_, or through the package manager in most Linux distributions. To run simulations in parallel using MPI, it is recommended to diff --git a/examples/xml/custom_source/CMakeLists.txt b/examples/xml/custom_source/CMakeLists.txt new file mode 100644 index 000000000..949817694 --- /dev/null +++ b/examples/xml/custom_source/CMakeLists.txt @@ -0,0 +1,8 @@ +cmake_minimum_required(VERSION 3.3 FATAL_ERROR) +project(openmc_sources CXX) +add_library(source SHARED source_ring.cpp) +find_package(OpenMC REQUIRED) +if (OpenMC_FOUND) + message(STATUS "Found OpenMC: ${OpenMC_DIR}") +endif() +target_link_libraries(source OpenMC::libopenmc) diff --git a/examples/xml/custom_source/geometry.xml b/examples/xml/custom_source/geometry.xml new file mode 100644 index 000000000..b30884f8c --- /dev/null +++ b/examples/xml/custom_source/geometry.xml @@ -0,0 +1,15 @@ + + + + + + + + + + + + + + + diff --git a/examples/xml/custom_source/materials.xml b/examples/xml/custom_source/materials.xml new file mode 100644 index 000000000..606c676df --- /dev/null +++ b/examples/xml/custom_source/materials.xml @@ -0,0 +1,16 @@ + + + + + + + + + + + + + + + + diff --git a/examples/xml/custom_source/settings.xml b/examples/xml/custom_source/settings.xml new file mode 100644 index 000000000..f935ed685 --- /dev/null +++ b/examples/xml/custom_source/settings.xml @@ -0,0 +1,14 @@ + + + + fixed source + 10 + 0 + 100000 + + + + build/libsource.so + + + diff --git a/examples/xml/custom_source/source_ring.cpp b/examples/xml/custom_source/source_ring.cpp new file mode 100644 index 000000000..d68122dd7 --- /dev/null +++ b/examples/xml/custom_source/source_ring.cpp @@ -0,0 +1,26 @@ +#include // for M_PI + +#include "openmc/random_lcg.h" +#include "openmc/source.h" +#include "openmc/particle.h" + +// you must have external C linkage here otherwise +// dlopen will not find the file +extern "C" openmc::Particle::Bank sample_source(uint64_t* seed) +{ + openmc::Particle::Bank particle; + // wgt + particle.particle = openmc::Particle::Type::neutron; + particle.wgt = 1.0; + // position + double angle = 2. * M_PI * openmc::prn(seed); + double radius = 3.0; + particle.r.x = radius * std::cos(angle); + particle.r.y = radius * std::sin(angle); + particle.r.z = 0.0; + // angle + particle.u = {1.0, 0.0, 0.0}; + particle.E = 14.08e6; + particle.delayed_group = 0; + return particle; +} diff --git a/examples/xml/custom_source/tallies.xml b/examples/xml/custom_source/tallies.xml new file mode 100644 index 000000000..7f6f29926 --- /dev/null +++ b/examples/xml/custom_source/tallies.xml @@ -0,0 +1,17 @@ + + + + + 100 + + + + 0 20.0e6 + + + + 1 2 + flux + + + diff --git a/include/openmc/settings.h b/include/openmc/settings.h index 2b4203646..494cfd1a2 100644 --- a/include/openmc/settings.h +++ b/include/openmc/settings.h @@ -60,6 +60,7 @@ extern std::string path_input; //!< directory where main .xml files r extern std::string path_output; //!< directory where output files are written extern std::string path_particle_restart; //!< path to a particle restart file extern std::string path_source; +extern std::string path_source_library; //!< path to the source shared object extern std::string path_sourcepoint; //!< path to a source file extern "C" std::string path_statepoint; //!< path to a statepoint file diff --git a/include/openmc/source.h b/include/openmc/source.h index a177995ea..4db17f241 100644 --- a/include/openmc/source.h +++ b/include/openmc/source.h @@ -68,6 +68,9 @@ Particle::Bank sample_external_source(uint64_t* seed); //! Fill source bank at end of generation for fixed source simulations void fill_source_bank_fixedsource(); +//! Fill source bank at the end of a generation for dlopen based source simulation +void fill_source_bank_custom_source(); + void free_memory_source(); } // namespace openmc diff --git a/openmc/data/photon.py b/openmc/data/photon.py index 32ed3010a..9cc697b78 100644 --- a/openmc/data/photon.py +++ b/openmc/data/photon.py @@ -784,8 +784,9 @@ class IncidentPhoton(EqualityMixin): sub_group = shell_group.create_group(key) # Write atomic relaxation - if key in self.atomic_relaxation.subshells: - self.atomic_relaxation.to_hdf5(sub_group, key) + if self.atomic_relaxation is not None: + if key in self.atomic_relaxation.subshells: + self.atomic_relaxation.to_hdf5(sub_group, key) else: continue diff --git a/openmc/deplete/nuclide.py b/openmc/deplete/nuclide.py index 7e352fbf0..0aae5adb7 100644 --- a/openmc/deplete/nuclide.py +++ b/openmc/deplete/nuclide.py @@ -533,6 +533,7 @@ class FissionYield(Mapping): return zip(self.products, self.yields) def __add__(self, other): + """Add one set of fission yields to this set, return new yields""" if not isinstance(other, FissionYield): return NotImplemented new = FissionYield(self.products, self.yields.copy()) @@ -550,12 +551,14 @@ class FissionYield(Mapping): return self + other def __imul__(self, scalar): + """Scale these fission yields by a real scalar""" if not isinstance(scalar, Real): return NotImplemented self.yields *= scalar return self def __mul__(self, scalar): + """Return a new set of yields scaled by a real scalar""" if not isinstance(scalar, Real): return NotImplemented new = FissionYield(self.products, self.yields.copy()) @@ -568,3 +571,8 @@ class FissionYield(Mapping): def __repr__(self): return "<{} containing {} products and yields>".format( self.__class__.__name__, len(self)) + + # Avoid greedy numpy operations like np.float64 * fission_yield + # converting this to an array on the fly. Force __rmul__ and + # __radd__. See issue #1492 + __array_ufunc__ = None diff --git a/openmc/model/funcs.py b/openmc/model/funcs.py index 3ee505240..9999c383f 100644 --- a/openmc/model/funcs.py +++ b/openmc/model/funcs.py @@ -373,52 +373,7 @@ def get_hexagonal_prism(*args, **kwargs): return hexagonal_prism(*args, **kwargs) -def cylinder_from_points(p1, p2, r, **kwargs): - """Return cylinder defined by two points passing through its center. - - Parameters - ---------- - p1, p2 : 3-tuples - Coordinates of two points that pass through the center of the cylinder - r : float - Radius of the cylinder - kwargs : dict - Keyword arguments passed to the :class:`openmc.Quadric` constructor - - Returns - ------- - openmc.Quadric - Quadric surface representing the cylinder. - - """ - # Get x, y, z coordinates of two points - x1, y1, z1 = p1 - x2, y2, z2 = p2 - - # Define intermediate terms - dx = x2 - x1 - dy = y2 - y1 - dz = z2 - z1 - cx = y1*z2 - y2*z1 - cy = x2*z1 - x1*z2 - cz = x1*y2 - x2*y1 - - # Given p=(x,y,z), p1=(x1, y1, z1), p2=(x2, y2, z2), the equation for the - # cylinder can be derived as r = |(p - p1) ⨯ (p - p2)| / |p2 - p1|. - # Expanding out all terms and grouping according to what Quadric expects - # gives the following coefficients. - kwargs['a'] = dy*dy + dz*dz - kwargs['b'] = dx*dx + dz*dz - kwargs['c'] = dx*dx + dy*dy - kwargs['d'] = -2*dx*dy - kwargs['e'] = -2*dy*dz - kwargs['f'] = -2*dx*dz - kwargs['g'] = 2*(cy*dz - cz*dy) - kwargs['h'] = 2*(cz*dx - cx*dz) - kwargs['j'] = 2*(cx*dy - cy*dx) - kwargs['k'] = cx*cx + cy*cy + cz*cz - (dx*dx + dy*dy + dz*dz)*r*r - - return Quadric(**kwargs) +cylinder_from_points = Cylinder.from_points def subdivide(surfaces): diff --git a/openmc/source.py b/openmc/source.py index 88c2f8611..7b709321f 100644 --- a/openmc/source.py +++ b/openmc/source.py @@ -8,20 +8,22 @@ from openmc.stats.multivariate import UnitSphere, Spatial import openmc.checkvalue as cv -class Source(object): +class Source: """Distribution of phase space coordinates for source sites. Parameters ---------- - space : openmc.stats.Spatial, optional + space : openmc.stats.Spatial Spatial distribution of source sites - angle : openmc.stats.UnitSphere, optional + angle : openmc.stats.UnitSphere Angular distribution of source sites - energy : openmc.stats.Univariate, optional + energy : openmc.stats.Univariate Energy distribution of source sites - filename : str, optional + filename : str Source file from which sites should be sampled - strength : Real + library : str + Path to a custom source library + strength : float Strength of the source particle : {'neutron', 'photon'} Source particle type @@ -36,7 +38,9 @@ class Source(object): Energy distribution of source sites file : str or None Source file from which sites should be sampled - strength : Real + library : str or None + Path to a custom source library + strength : float Strength of the source particle : {'neutron', 'photon'} Source particle type @@ -44,11 +48,12 @@ class Source(object): """ def __init__(self, space=None, angle=None, energy=None, filename=None, - strength=1.0, particle='neutron'): + library=None, strength=1.0, particle='neutron'): self._space = None self._angle = None self._energy = None self._file = None + self._library = None if space is not None: self.space = space @@ -58,6 +63,8 @@ class Source(object): self.energy = energy if filename is not None: self.file = filename + if library is not None: + self.library = library self.strength = strength self.particle = particle @@ -65,6 +72,10 @@ class Source(object): def file(self): return self._file + @property + def library(self): + return self._library + @property def space(self): return self._space @@ -90,6 +101,11 @@ class Source(object): cv.check_type('source file', filename, str) self._file = filename + @library.setter + def library(self, library_name): + cv.check_type('library', library_name, str) + self._library = library_name + @space.setter def space(self, space): cv.check_type('spatial distribution', space, Spatial) @@ -131,6 +147,8 @@ class Source(object): element.set("particle", self.particle) if self.file is not None: element.set("file", self.file) + if self.library is not None: + element.set("library", self.library) if self.space is not None: element.append(self.space.to_xml_element()) if self.angle is not None: @@ -168,6 +186,10 @@ class Source(object): if filename is not None: source.file = filename + library = get_text(elem, 'library') + if library is not None: + source.library = library + space = elem.find('space') if space is not None: source.space = Spatial.from_xml_element(space) diff --git a/openmc/surface.py b/openmc/surface.py index 55e3f459d..567968b2c 100644 --- a/openmc/surface.py +++ b/openmc/surface.py @@ -1,15 +1,16 @@ from abc import ABCMeta, abstractmethod from collections import OrderedDict from copy import deepcopy -from numbers import Real, Integral +from numbers import Real from xml.etree import ElementTree as ET -from warnings import warn +from warnings import warn, catch_warnings, simplefilter +import math import numpy as np from openmc.checkvalue import check_type, check_value from openmc.region import Region, Intersection, Union -from openmc.mixin import IDManagerMixin +from openmc.mixin import IDManagerMixin, IDWarning _BOUNDARY_TYPES = ['transmission', 'vacuum', 'reflective', 'periodic', 'white'] @@ -19,6 +20,11 @@ _WARNING_UPPER = """\ will not accept the capitalized version.\ """ +_WARNING_KWARGS = """\ +"{}(...) accepts keyword arguments only for '{}'. Future versions of OpenMC \ +will not accept positional parameters for superclass arguments.\ +""" + class Surface(IDManagerMixin, metaclass=ABCMeta): """An implicit surface with an associated boundary condition. @@ -60,6 +66,7 @@ class Surface(IDManagerMixin, metaclass=ABCMeta): next_id = 1 used_ids = set() + _atol = 1.e-12 def __init__(self, surface_id=None, boundary_type='transmission', name=''): self.id = surface_id @@ -179,13 +186,85 @@ class Surface(IDManagerMixin, metaclass=ABCMeta): return memo[self] - @abstractmethod - def evaluate(self, point): - pass + def normalize(self, coeffs=None): + """Normalize coefficients by first nonzero value + + Parameters + ---------- + coeffs : tuple, optional + Tuple of surface coefficients to normalize. Defaults to None. If no + coefficients are supplied then the coefficients will be taken from + the current Surface. + + Returns + ------- + tuple of normalized coefficients + + """ + if coeffs is None: + coeffs = self._get_base_coeffs() + coeffs = np.asarray(coeffs) + nonzeros = ~np.isclose(coeffs, 0., rtol=0., atol=self._atol) + norm_factor = np.abs(coeffs[nonzeros][0]) + return tuple([c/norm_factor for c in coeffs]) + + def is_equal(self, other): + """Determine if this Surface is equivalent to another + + Parameters + ---------- + other : instance of openmc.Surface + Instance of openmc.Surface that should be compared to the current + surface + + """ + coeffs1 = self.normalize(self._get_base_coeffs()) + coeffs2 = self.normalize(other._get_base_coeffs()) + + return np.allclose(coeffs1, coeffs2, rtol=0., atol=self._atol) @abstractmethod - def translate(self, vector): - pass + def _get_base_coeffs(self): + """Return polynomial coefficients representing the implicit surface + equation. + + """ + + @abstractmethod + def evaluate(self, point): + """Evaluate the surface equation at a given point. + + Parameters + ---------- + point : 3-tuple of float + The Cartesian coordinates, :math:`(x',y',z')`, at which the surface + equation should be evaluated. + + Returns + ------- + float + Evaluation of the surface polynomial at point :math:`(x',y',z')` + + """ + + @abstractmethod + def translate(self, vector, inplace=False): + """Translate surface in given direction + + Parameters + ---------- + vector : iterable of float + Direction in which surface should be translated + inplace : boolean + Whether or not to return a new instance of this Surface or to + modify the coefficients of this Surface. Defaults to False + + Returns + ------- + instance of openmc.Surface + Translated surface + + """ def to_xml_element(self): """Return XML representation of the surface @@ -228,21 +307,7 @@ class Surface(IDManagerMixin, metaclass=ABCMeta): # Determine appropriate class surf_type = elem.get('type') - surface_classes = { - 'plane': Plane, - 'x-plane': XPlane, - 'y-plane': YPlane, - 'z-plane': ZPlane, - 'x-cylinder': XCylinder, - 'y-cylinder': YCylinder, - 'z-cylinder': ZCylinder, - 'sphere': Sphere, - 'x-cone': XCone, - 'y-cone': YCone, - 'z-cone': ZCone, - 'quadric': Quadric, - } - cls = surface_classes[surf_type] + cls = _SURFACE_CLASSES[surf_type] # Determine ID, boundary type, coefficients kwargs = {} @@ -268,62 +333,155 @@ class Surface(IDManagerMixin, metaclass=ABCMeta): Instance of surface subclass """ + surface_id = int(group.name.split('/')[-1].lstrip('surface ')) name = group['name'][()].decode() if 'name' in group else '' surf_type = group['type'][()].decode() bc = group['boundary_type'][()].decode() coeffs = group['coefficients'][...] + kwargs = {'boundary_type': bc, 'name': name, 'surface_id': surface_id} - # Create the Surface based on its type - if surf_type == 'x-plane': - x0 = coeffs[0] - surface = XPlane(x0, bc, name, surface_id) + cls = _SURFACE_CLASSES[surf_type] - elif surf_type == 'y-plane': - y0 = coeffs[0] - surface = YPlane(y0, bc, name, surface_id) - - elif surf_type == 'z-plane': - z0 = coeffs[0] - surface = ZPlane(z0, bc, name, surface_id) - - elif surf_type == 'plane': - A, B, C, D = coeffs - surface = Plane(A, B, C, D, bc, name, surface_id) - - elif surf_type == 'x-cylinder': - y0, z0, r = coeffs - surface = XCylinder(y0, z0, r, bc, name, surface_id) - - elif surf_type == 'y-cylinder': - x0, z0, r = coeffs - surface = YCylinder(x0, z0, r, bc, name, surface_id) - - elif surf_type == 'z-cylinder': - x0, y0, r = coeffs - surface = ZCylinder(x0, y0, r, bc, name, surface_id) - - elif surf_type == 'sphere': - x0, y0, z0, r = coeffs - surface = Sphere(x0, y0, z0, r, bc, name, surface_id) - - elif surf_type in ['x-cone', 'y-cone', 'z-cone']: - x0, y0, z0, r2 = coeffs - if surf_type == 'x-cone': - surface = XCone(x0, y0, z0, r2, bc, name, surface_id) - elif surf_type == 'y-cone': - surface = YCone(x0, y0, z0, r2, bc, name, surface_id) - elif surf_type == 'z-cone': - surface = ZCone(x0, y0, z0, r2, bc, name, surface_id) - - elif surf_type == 'quadric': - a, b, c, d, e, f, g, h, j, k = coeffs - surface = Quadric(a, b, c, d, e, f, g, h, j, k, bc, name, surface_id) - - return surface + return cls(*coeffs, **kwargs) -class Plane(Surface): +class PlaneMixin(metaclass=ABCMeta): + """A Plane mixin class for all operations on order 1 surfaces""" + def __init__(self, **kwargs): + super().__init__(**kwargs) + self._periodic_surface = None + + @property + def periodic_surface(self): + return self._periodic_surface + + @periodic_surface.setter + def periodic_surface(self, periodic_surface): + check_type('periodic surface', periodic_surface, Plane) + self._periodic_surface = periodic_surface + periodic_surface._periodic_surface = self + + def _get_base_coeffs(self): + return (self.a, self.b, self.c, self.d) + + def _get_normal(self): + a, b, c = self._get_base_coeffs()[:3] + return np.array((a, b, c)) / math.sqrt(a*a + b*b + c*c) + + def bounding_box(self, side): + """Determine an axis-aligned bounding box. + + An axis-aligned bounding box for Plane half-spaces is represented by + its lower-left and upper-right coordinates. If the half-space is + unbounded in a particular direction, numpy.inf is used to represent + infinity. + + Parameters + ---------- + side : {'+', '-'} + Indicates the negative or positive half-space + + Returns + ------- + numpy.ndarray + Lower-left coordinates of the axis-aligned bounding box for the + desired half-space + numpy.ndarray + Upper-right coordinates of the axis-aligned bounding box for the + desired half-space + + """ + # Compute the bounding box based on the normal vector to the plane + nhat = self._get_normal() + ll = np.array([-np.inf, -np.inf, -np.inf]) + ur = np.array([np.inf, np.inf, np.inf]) + # If the plane is axis aligned, find the proper bounding box + if np.any(np.isclose(np.abs(nhat), 1., rtol=0., atol=self._atol)): + sign = nhat.sum() + a, b, c, d = self._get_base_coeffs() + vals = [d/val if not np.isclose(val, 0., rtol=0., atol=self._atol) + else np.nan for val in (a, b, c)] + if side == '-': + if sign > 0: + ur = np.array([v if not np.isnan(v) else np.inf for v in vals]) + else: + ll = np.array([v if not np.isnan(v) else -np.inf for v in vals]) + elif side == '+': + if sign > 0: + ll = np.array([v if not np.isnan(v) else -np.inf for v in vals]) + else: + ur = np.array([v if not np.isnan(v) else np.inf for v in vals]) + + return (ll, ur) + + def evaluate(self, point): + """Evaluate the surface equation at a given point. + + Parameters + ---------- + point : 3-tuple of float + The Cartesian coordinates, :math:`(x',y',z')`, at which the surface + equation should be evaluated. + + Returns + ------- + float + :math:`Ax' + By' + Cz' - D` + + """ + + x, y, z = point + a, b, c, d = self._get_base_coeffs() + return a*x + b*y + c*z - d + + def translate(self, vector, inplace=False): + """Translate surface in given direction + + Parameters + ---------- + vector : iterable of float + Direction in which surface should be translated + inplace : boolean + Whether or not to return a new instance of a Plane or to modify the + coefficients of this plane. Defaults to False + + Returns + ------- + openmc.Plane + Translated surface + + """ + vx, vy, vz = vector + a, b, c, d = self._get_base_coeffs() + d = d + a*vx + b*vy + c*vz + + surf = self if inplace else self.clone() + + setattr(surf, surf._coeff_keys[-1], d) + + return surf + + def to_xml_element(self): + """Return XML representation of the surface + + Returns + ------- + element : xml.etree.ElementTree.Element + XML element containing source data + + """ + element = super().to_xml_element() + + # Add periodic surface pair information + if self.boundary_type == 'periodic': + if self.periodic_surface is not None: + element.set("periodic_surface_id", + str(self.periodic_surface.id)) + return element + + +class Plane(PlaneMixin, Surface): """An arbitrary plane of the form :math:`Ax + By + Cz = D`. Parameters @@ -376,19 +534,40 @@ class Plane(Surface): _type = 'plane' _coeff_keys = ('a', 'b', 'c', 'd') - def __init__(self, a=1., b=0., c=0., d=0., boundary_type='transmission', - name='', surface_id=None, **kwargs): - super().__init__(surface_id, boundary_type, name=name) - self._periodic_surface = None - self.a = a - self.b = b - self.c = c - self.d = d - for k, v in kwargs.items(): - if k in 'ABCD': - warn(_WARNING_UPPER.format(type(self).__name__, k.lower(), k), - FutureWarning) - setattr(self, k.lower(), v) + def __init__(self, a=1., b=0., c=0., d=0., *args, **kwargs): + # *args should ultimately be limited to a, b, c, d as specified in + # __init__, but to preserve the API it is allowed to accept Surface + # parameters for now, but will raise warnings if this is done. + argtup = ('a', 'b', 'c', 'd', 'boundary_type', 'name', 'surface_id') + kwargs.update(dict(zip(argtup, args))) + + # Warn if Surface parameters are passed by position, not by keyword + superkwargs = {} + for k in ('boundary_type', 'name', 'surface_id'): + val = kwargs.get(k, None) + if val is not None: + superkwargs[k] = val + warn(_WARNING_KWARGS.format(type(self), k), + FutureWarning) + + super().__init__(**superkwargs) + + for key, val in zip(self._coeff_keys, (a, b, c, d)): + setattr(self, key, val) + + # Warn if capital letter arguments are passed + for k in 'ABCD': + val = kwargs.pop(k, None) + if val is not None: + warn(_WARNING_UPPER.format(type(self), k.lower(), k), + FutureWarning) + setattr(self, k.lower(), val) + + @classmethod + def __subclasshook__(cls, c): + if cls is Plane and c in (XPlane, YPlane, ZPlane): + return True + return NotImplemented @property def a(self): @@ -406,10 +585,6 @@ class Plane(Surface): def d(self): return self.coefficients['d'] - @property - def periodic_surface(self): - return self._periodic_surface - @a.setter def a(self, a): check_type('A coefficient', a, Real) @@ -430,69 +605,6 @@ class Plane(Surface): check_type('D coefficient', d, Real) self._coefficients['d'] = d - @periodic_surface.setter - def periodic_surface(self, periodic_surface): - check_type('periodic surface', periodic_surface, Plane) - self._periodic_surface = periodic_surface - periodic_surface._periodic_surface = self - - def evaluate(self, point): - """Evaluate the surface equation at a given point. - - Parameters - ---------- - point : 3-tuple of float - The Cartesian coordinates, :math:`(x',y',z')`, at which the surface - equation should be evaluated. - - Returns - ------- - float - :math:`Ax' + By' + Cz' - D` - - """ - - x, y, z = point - return self.a*x + self.b*y + self.c*z - self.d - - def translate(self, vector): - """Translate surface in given direction - - Parameters - ---------- - vector : iterable of float - Direction in which surface should be translated - - Returns - ------- - openmc.Plane - Translated surface - - """ - vx, vy, vz = vector - d = self.d + self.a*vx + self.b*vy + self.c*vz - if d == self.d: - return self - else: - return type(self)(a=self.a, b=self.b, c=self.c, d=d) - - def to_xml_element(self): - """Return XML representation of the surface - - Returns - ------- - element : xml.etree.ElementTree.Element - XML element containing source data - - """ - element = super().to_xml_element() - - # Add periodic surface pair information - if self.boundary_type == 'periodic': - if self.periodic_surface is not None: - element.set("periodic_surface_id", str(self.periodic_surface.id)) - return element - @classmethod def from_points(cls, p1, p2, p3, **kwargs): """Return a plane given three points that pass through it. @@ -526,7 +638,7 @@ class Plane(Surface): return cls(a=a, b=b, c=c, d=d, **kwargs) -class XPlane(Plane): +class XPlane(PlaneMixin, Surface): """A plane perpendicular to the x axis of the form :math:`x - x_0 = 0` Parameters @@ -568,90 +680,47 @@ class XPlane(Plane): _type = 'x-plane' _coeff_keys = ('x0',) - def __init__(self, x0=0., boundary_type='transmission', - name='', surface_id=None): - super().__init__(surface_id=surface_id, boundary_type=boundary_type, name=name) + def __init__(self, x0=0., *args, **kwargs): + # work around for accepting Surface kwargs as positional parameters + # until they are deprecated + argsdict = dict(zip(('boundary_type', 'name', 'surface_id'), args)) + for k in argsdict: + warn(_WARNING_KWARGS.format(type(self).__name__, k), FutureWarning) + kwargs.update(argsdict) + + super().__init__(**kwargs) self.x0 = x0 @property def x0(self): return self.coefficients['x0'] + @property + def a(self): + return 1. + + @property + def b(self): + return 0. + + @property + def c(self): + return 0. + + @property + def d(self): + return self.x0 + @x0.setter def x0(self, x0): check_type('x0 coefficient', x0, Real) self._coefficients['x0'] = x0 - def bounding_box(self, side): - """Determine an axis-aligned bounding box. - - An axis-aligned bounding box for surface half-spaces is represented by - its lower-left and upper-right coordinates. For the x-plane surface, the - half-spaces are unbounded in their y- and z- directions. To represent - infinity, numpy.inf is used. - - Parameters - ---------- - side : {'+', '-'} - Indicates the negative or positive half-space - - Returns - ------- - numpy.ndarray - Lower-left coordinates of the axis-aligned bounding box for the - desired half-space - numpy.ndarray - Upper-right coordinates of the axis-aligned bounding box for the - desired half-space - - """ - - if side == '-': - return (np.array([-np.inf, -np.inf, -np.inf]), - np.array([self.x0, np.inf, np.inf])) - elif side == '+': - return (np.array([self.x0, -np.inf, -np.inf]), - np.array([np.inf, np.inf, np.inf])) - def evaluate(self, point): - """Evaluate the surface equation at a given point. - - Parameters - ---------- - point : 3-tuple of float - The Cartesian coordinates, :math:`(x',y',z')`, at which the surface - equation should be evaluated. - - Returns - ------- - float - :math:`x' - x_0` - - """ return point[0] - self.x0 - def translate(self, vector): - """Translate surface in given direction - Parameters - ---------- - vector : iterable of float - Direction in which surface should be translated - - Returns - ------- - openmc.XPlane - Translated surface - - """ - vx = vector[0] - if vx == 0: - return self - else: - return type(self)(x0=self.x0 + vx) - - -class YPlane(Plane): +class YPlane(PlaneMixin, Surface): """A plane perpendicular to the y axis of the form :math:`y - y_0 = 0` Parameters @@ -693,106 +762,63 @@ class YPlane(Plane): _type = 'y-plane' _coeff_keys = ('y0',) - def __init__(self, y0=0., boundary_type='transmission', - name='', surface_id=None): - super().__init__(surface_id=surface_id, boundary_type=boundary_type, name=name) + def __init__(self, y0=0., *args, **kwargs): + # work around for accepting Surface kwargs as positional parameters + # until they are deprecated + argsdict = dict(zip(('boundary_type', 'name', 'surface_id'), args)) + for k in argsdict: + warn(_WARNING_KWARGS.format(type(self).__name__, k), FutureWarning) + kwargs.update(argsdict) + + super().__init__(**kwargs) self.y0 = y0 @property def y0(self): return self.coefficients['y0'] + @property + def a(self): + return 0. + + @property + def b(self): + return 1. + + @property + def c(self): + return 0. + + @property + def d(self): + return self.y0 + @y0.setter def y0(self, y0): check_type('y0 coefficient', y0, Real) self._coefficients['y0'] = y0 - def bounding_box(self, side): - """Determine an axis-aligned bounding box. - - An axis-aligned bounding box for surface half-spaces is represented by - its lower-left and upper-right coordinates. For the y-plane surface, the - half-spaces are unbounded in their x- and z- directions. To represent - infinity, numpy.inf is used. - - Parameters - ---------- - side : {'+', '-'} - Indicates the negative or positive half-space - - Returns - ------- - numpy.ndarray - Lower-left coordinates of the axis-aligned bounding box for the - desired half-space - numpy.ndarray - Upper-right coordinates of the axis-aligned bounding box for the - desired half-space - - """ - - if side == '-': - return (np.array([-np.inf, -np.inf, -np.inf]), - np.array([np.inf, self.y0, np.inf])) - elif side == '+': - return (np.array([-np.inf, self.y0, -np.inf]), - np.array([np.inf, np.inf, np.inf])) - def evaluate(self, point): - """Evaluate the surface equation at a given point. - - Parameters - ---------- - point : 3-tuple of float - The Cartesian coordinates, :math:`(x',y',z')`, at which the surface - equation should be evaluated. - - Returns - ------- - float - :math:`y' - y_0` - - """ return point[1] - self.y0 - def translate(self, vector): - """Translate surface in given direction - Parameters - ---------- - vector : iterable of float - Direction in which surface should be translated - - Returns - ------- - openmc.YPlane - Translated surface - - """ - vy = vector[1] - if vy == 0.0: - return self - else: - return type(self)(y0=self.y0 + vy) - - -class ZPlane(Plane): +class ZPlane(PlaneMixin, Surface): """A plane perpendicular to the z axis of the form :math:`z - z_0 = 0` Parameters ---------- - surface_id : int, optional - Unique identifier for the surface. If not specified, an identifier will - automatically be assigned. + z0 : float, optional + Location of the plane. Defaults to 0. boundary_type : {'transmission, 'vacuum', 'reflective', 'periodic', 'white'}, optional Boundary condition that defines the behavior for particles hitting the surface. Defaults to transmissive boundary condition where particles freely pass through the surface. Only axis-aligned periodicity is supported, i.e., x-planes can only be paired with x-planes. - z0 : float, optional - Location of the plane. Defaults to 0. name : str, optional Name of the plane. If not specified, the name will be the empty string. + surface_id : int, optional + Unique identifier for the surface. If not specified, an identifier will + automatically be assigned. Attributes ---------- @@ -818,50 +844,99 @@ class ZPlane(Plane): _type = 'z-plane' _coeff_keys = ('z0',) - def __init__(self, z0=0., boundary_type='transmission', - name='', surface_id=None): - super().__init__(surface_id=surface_id, boundary_type=boundary_type, name=name) + def __init__(self, z0=0., *args, **kwargs): + # work around for accepting Surface kwargs as positional parameters + # until they are deprecated + argsdict = dict(zip(('boundary_type', 'name', 'surface_id'), args)) + for k in argsdict: + warn(_WARNING_KWARGS.format(type(self).__name__, k), FutureWarning) + kwargs.update(argsdict) + + super().__init__(**kwargs) self.z0 = z0 @property def z0(self): return self.coefficients['z0'] + @property + def a(self): + return 0. + + @property + def b(self): + return 0. + + @property + def c(self): + return 1. + + @property + def d(self): + return self.z0 + @z0.setter def z0(self, z0): check_type('z0 coefficient', z0, Real) self._coefficients['z0'] = z0 - def bounding_box(self, side): - """Determine an axis-aligned bounding box. + def evaluate(self, point): + return point[2] - self.z0 - An axis-aligned bounding box for surface half-spaces is represented by - its lower-left and upper-right coordinates. For the z-plane surface, the - half-spaces are unbounded in their x- and y- directions. To represent - infinity, numpy.inf is used. + +class QuadricMixin(metaclass=ABCMeta): + """A Mixin class implementing common functionality for quadric surfaces""" + + @property + def _origin(self): + return np.array((self.x0, self.y0, self.z0)) + + @property + def _axis(self): + axis = np.array((self.dx, self.dy, self.dz)) + return axis / np.linalg.norm(axis) + + def get_Abc(self, coeffs=None): + """Compute matrix, vector, and scalar coefficients for this surface or + for a specified set of coefficients. Parameters ---------- - side : {'+', '-'} - Indicates the negative or positive half-space + coeffs : tuple, optional + Tuple of coefficients from which to compute the quadric elements. + If none are supplied the coefficients of this surface will be used. + """ + if coeffs is None: + a, b, c, d, e, f, g, h, j, k = self._get_base_coeffs() + else: + a, b, c, d, e, f, g, h, j, k = coeffs + + A = np.array([[a, d/2, f/2], [d/2, b, e/2], [f/2, e/2, c]]) + bvec = np.array([g, h, j]) + + return A, bvec, k + + def eigh(self, coeffs=None): + """Wrapper method for returning eigenvalues and eigenvectors of this + quadric surface which is used for transformations. + + Parameters + ---------- + coeffs : tuple, optional + Tuple of coefficients from which to compute the quadric elements. + If none are supplied the coefficients of this surface will be used. Returns ------- - numpy.ndarray - Lower-left coordinates of the axis-aligned bounding box for the - desired half-space - numpy.ndarray - Upper-right coordinates of the axis-aligned bounding box for the - desired half-space + w, v : tuple of numpy arrays with shapes (3,) and (3,3) respectively + Returns the eigenvalues and eigenvectors of the quadric matrix A + that represents the supplied coefficients. The vector w contains + the eigenvalues in ascending order and the matrix v contains the + eigenvectors such that v[:,i] is the eigenvector corresponding to + the eigenvalue w[i]. """ - - if side == '-': - return (np.array([-np.inf, -np.inf, -np.inf]), - np.array([np.inf, np.inf, self.z0])) - elif side == '+': - return (np.array([-np.inf, -np.inf, self.z0]), - np.array([np.inf, np.inf, np.inf])) + return np.linalg.eigh(self.get_Abc(coeffs=coeffs)[0]) def evaluate(self, point): """Evaluate the surface equation at a given point. @@ -875,37 +950,287 @@ class ZPlane(Plane): Returns ------- float - :math:`z' - z_0` + :math:`Ax'^2 + By'^2 + Cz'^2 + Dx'y' + Ey'z' + Fx'z' + Gx' + Hy' + + Jz' + K = 0` """ - return point[2] - self.z0 + x = np.asarray(point) + A, b, c = self.get_Abc() + return x.T @ A @ x + b.T @ x + c - def translate(self, vector): + def translate(self, vector, inplace=False): """Translate surface in given direction Parameters ---------- vector : iterable of float Direction in which surface should be translated + inplace : boolean + Whether to return a clone of the Surface or the Surface itself. + Defaults to False Returns ------- - openmc.ZPlane + openmc.Surface Translated surface """ - vz = vector[2] - if vz == 0.0: - return self + vector = np.asarray(vector) + + surf = self if inplace else self.clone() + + if hasattr(self, 'x0'): + for vi, xi in zip(vector, ('x0', 'y0', 'z0')): + val = getattr(surf, xi) + try: + setattr(surf, xi, val + vi) + except AttributeError: + # That attribute is read only i.e x0 for XCylinder + pass + else: - return type(self)(z0=self.z0 + vz) + A, bvec, cnst = self.get_Abc() + + g, h, j = bvec - 2*vector.T @ A + k = cnst + vector.T @ A @ vector - bvec.T @ vector + + for key, val in zip(('g', 'h', 'j', 'k'), (g, h, j, k)): + setattr(surf, key, val) + + return surf -class Cylinder(Surface): - """A cylinder whose length is parallel to the x-, y-, or z-axis. +class Cylinder(QuadricMixin, Surface): + """A cylinder with radius r, centered on the point (x0, y0, z0) with an + axis specified by the line through points (x0, y0, z0) and (x0+dx, y0+dy, + z0+dz) Parameters ---------- + x0 : float, optional + x-coordinate for the origin of the Cylinder. Defaults to 0 + y0 : float, optional + y-coordinate for the origin of the Cylinder. Defaults to 0 + z0 : float, optional + z-coordinate for the origin of the Cylinder. Defaults to 0 + r : float, optional + Radius of the cylinder. Defaults to 1. + dx : float, optional + x-component of the vector representing the axis of the cylinder. + Defaults to 0. + dy : float, optional + y-component of the vector representing the axis of the cylinder. + Defaults to 0. + dz : float, optional + z-component of the vector representing the axis of the cylinder. + Defaults to 1. + boundary_type : {'transmission, 'vacuum', 'reflective', 'white'}, optional + Boundary condition that defines the behavior for particles hitting the + surface. Defaults to transmissive boundary condition where particles + freely pass through the surface. + name : str, optional + Name of the cylinder. If not specified, the name will be the empty + string. + surface_id : int, optional + Unique identifier for the surface. If not specified, an identifier will + automatically be assigned. + + Attributes + ---------- + x0 : float + x-coordinate for the origin of the Cylinder + y0 : float + y-coordinate for the origin of the Cylinder + z0 : float + z-coordinate for the origin of the Cylinder + r : float + Radius of the cylinder + dx : float + x-component of the vector representing the axis of the cylinder + dy : float + y-component of the vector representing the axis of the cylinder + dz : float + z-component of the vector representing the axis of the cylinder + boundary_type : {'transmission, 'vacuum', 'reflective', 'white'} + Boundary condition that defines the behavior for particles hitting the + surface. + coefficients : dict + Dictionary of surface coefficients + id : int + Unique identifier for the surface + name : str + Name of the surface + type : str + Type of the surface + + """ + _type = 'cylinder' + _coeff_keys = ('x0', 'y0', 'z0', 'r', 'dx', 'dy', 'dz') + + def __init__(self, x0=0., y0=0., z0=0., r=1., dx=0., dy=0., dz=1., **kwargs): + super().__init__(**kwargs) + + for key, val in zip(self._coeff_keys, (x0, y0, z0, r, dx, dy, dz)): + setattr(self, key, val) + + @classmethod + def __subclasshook__(cls, c): + if cls is Cylinder and c in (XCylinder, YCylinder, ZCylinder): + return True + return NotImplemented + + @property + def x0(self): + return self.coefficients['x0'] + + @property + def y0(self): + return self.coefficients['y0'] + + @property + def z0(self): + return self.coefficients['z0'] + + @property + def r(self): + return self.coefficients['r'] + + @property + def dx(self): + return self.coefficients['dx'] + + @property + def dy(self): + return self.coefficients['dy'] + + @property + def dz(self): + return self.coefficients['dz'] + + @x0.setter + def x0(self, x0): + check_type('x0 coefficient', x0, Real) + self._coefficients['x0'] = x0 + + @y0.setter + def y0(self, y0): + check_type('y0 coefficient', y0, Real) + self._coefficients['y0'] = y0 + + @z0.setter + def z0(self, z0): + check_type('z0 coefficient', z0, Real) + self._coefficients['z0'] = z0 + + @r.setter + def r(self, r): + check_type('r coefficient', r, Real) + self._coefficients['r'] = r + + @dx.setter + def dx(self, dx): + check_type('dx coefficient', dx, Real) + self._coefficients['dx'] = dx + + @dy.setter + def dy(self, dy): + check_type('dy coefficient', dy, Real) + self._coefficients['dy'] = dy + + @dz.setter + def dz(self, dz): + check_type('dz coefficient', dz, Real) + self._coefficients['dz'] = dz + + def _get_base_coeffs(self): + # Get x, y, z coordinates of two points + x1, y1, z1 = self._origin + x2, y2, z2 = self._origin + self._axis + r = self.r + + # Define intermediate terms + dx = x2 - x1 + dy = y2 - y1 + dz = z2 - z1 + cx = y1*z2 - y2*z1 + cy = x2*z1 - x1*z2 + cz = x1*y2 - x2*y1 + + # Given p=(x,y,z), p1=(x1, y1, z1), p2=(x2, y2, z2), the equation + # for the cylinder can be derived as + # r = |(p - p1) ⨯ (p - p2)| / |p2 - p1|. + # Expanding out all terms and grouping according to what Quadric + # expects gives the following coefficients. + a = dy*dy + dz*dz + b = dx*dx + dz*dz + c = dx*dx + dy*dy + d = -2*dx*dy + e = -2*dy*dz + f = -2*dx*dz + g = 2*(cy*dz - cz*dy) + h = 2*(cz*dx - cx*dz) + j = 2*(cx*dy - cy*dx) + k = cx*cx + cy*cy + cz*cz - (dx*dx + dy*dy + dz*dz)*r*r + + return (a, b, c, d, e, f, g, h, j, k) + + @classmethod + def from_points(cls, p1, p2, r=1., **kwargs): + """Return a cylinder given points that define the axis and a radius. + + Parameters + ---------- + p1, p2 : 3-tuples + Points that pass through the plane, p1 will be used as (x0, y0, z0) + r : float, optional + Radius of the cylinder. Defaults to 1. + kwargs : dict + Keyword arguments passed to the :class:`Cylinder` constructor + + Returns + ------- + Cylinder + Cylinder that has an axis through the points p1 and p2, and a + radius r. + + """ + # Convert to numpy arrays + p1 = np.asarray(p1) + p2 = np.asarray(p2) + x0, y0, z0 = p1 + dx, dy, dz = p2 - p1 + + return cls(x0=x0, y0=y0, z0=z0, r=r, dx=dx, dy=dy, dz=dz, **kwargs) + + def to_xml_element(self): + """Return XML representation of the surface + + Returns + ------- + element : xml.etree.ElementTree.Element + XML element containing source data + + """ + # This method overrides Surface.to_xml_element to generate a Quadric + # since the C++ layer doesn't support Cylinders right now + with catch_warnings(): + simplefilter('ignore', IDWarning) + kwargs = {'boundary_type': self.boundary_type, 'name': self.name, + 'surface_id': self.id} + quad_rep = Quadric(*self._get_base_coeffs(), **kwargs) + return quad_rep.to_xml_element() + + +class XCylinder(QuadricMixin, Surface): + """An infinite cylinder whose length is parallel to the x-axis of the form + :math:`(y - y_0)^2 + (z - z_0)^2 = r^2`. + + Parameters + ---------- + y0 : float, optional + y-coordinate for the origin of the Cylinder. Defaults to 0 + z0 : float, optional + z-coordinate for the origin of the Cylinder. Defaults to 0 r : float, optional Radius of the cylinder. Defaults to 1. boundary_type : {'transmission, 'vacuum', 'reflective', 'white'}, optional @@ -921,6 +1246,10 @@ class Cylinder(Surface): Attributes ---------- + y0 : float + y-coordinate for the origin of the Cylinder + z0 : float + z-coordinate for the origin of the Cylinder r : float Radius of the cylinder boundary_type : {'transmission, 'vacuum', 'reflective', 'white'} @@ -935,76 +1264,21 @@ class Cylinder(Surface): type : str Type of the surface - """ - def __init__(self, r=1., boundary_type='transmission', - name='', surface_id=None): - super().__init__(surface_id, boundary_type, name=name) - self.r = r - - @property - def r(self): - return self.coefficients['r'] - - @r.setter - def r(self, r): - check_type('r coefficient', r, Real) - self._coefficients['r'] = r - - -class XCylinder(Cylinder): - """An infinite cylinder whose length is parallel to the x-axis of the form - :math:`(y - y_0)^2 + (z - z_0)^2 = r^2`. - - Parameters - ---------- - y0 : float, optional - y-coordinate of the center of the cylinder. Defaults to 0. - z0 : float, optional - z-coordinate of the center of the cylinder. Defaults to 0. - r : float, optional - Radius of the cylinder. Defaults to 0. - boundary_type : {'transmission, 'vacuum', 'reflective', 'white'}, optional - Boundary condition that defines the behavior for particles hitting the - surface. Defaults to transmissive boundary condition where particles - freely pass through the surface. - name : str, optional - Name of the cylinder. If not specified, the name will be the empty - string. - surface_id : int, optional - Unique identifier for the surface. If not specified, an identifier will - automatically be assigned. - - Attributes - ---------- - y0 : float - y-coordinate of the center of the cylinder - z0 : float - z-coordinate of the center of the cylinder - boundary_type : {'transmission, 'vacuum', 'reflective', 'white'} - Boundary condition that defines the behavior for particles hitting the - surface. - coefficients : dict - Dictionary of surface coefficients - id : int - Unique identifier for the surface - name : str - Name of the surface - type : str - Type of the surface - """ _type = 'x-cylinder' _coeff_keys = ('y0', 'z0', 'r') - def __init__(self, y0=0., z0=0., r=1., boundary_type='transmission', - name='', surface_id=None, *, R=None): + def __init__(self, y0=0., z0=0., r=1., **kwargs): + R = kwargs.pop('R', None) if R is not None: - warn(_WARNING_UPPER.format(type(self).__name__, 'r', 'R'), FutureWarning) + warn(_WARNING_UPPER.format(type(self).__name__, 'r', 'R'), + FutureWarning) r = R - super().__init__(r, boundary_type, name, surface_id) - self.y0 = y0 - self.z0 = z0 + super().__init__(**kwargs) + + for key, val in zip(self._coeff_keys, (y0, z0, r)): + setattr(self, key, val) @property def y0(self): @@ -1014,6 +1288,26 @@ class XCylinder(Cylinder): def z0(self): return self.coefficients['z0'] + @property + def r(self): + return self.coefficients['r'] + + @property + def x0(self): + return 0. + + @property + def dx(self): + return 1. + + @property + def dy(self): + return 0. + + @property + def dz(self): + return 0. + @y0.setter def y0(self, y0): check_type('y0 coefficient', y0, Real) @@ -1024,31 +1318,21 @@ class XCylinder(Cylinder): check_type('z0 coefficient', z0, Real) self._coefficients['z0'] = z0 + @r.setter + def r(self, r): + check_type('r coefficient', r, Real) + self._coefficients['r'] = r + + def _get_base_coeffs(self): + y0, z0, r = self.y0, self.z0, self.r + + a = d = e = f = g = 0. + b = c = 1. + h, j, k = -2*y0, -2*z0, y0*y0 + z0*z0 - r*r + + return (a, b, c, d, e, f, g, h, j, k) + def bounding_box(self, side): - """Determine an axis-aligned bounding box. - - An axis-aligned bounding box for surface half-spaces is represented by - its lower-left and upper-right coordinates. For the x-cylinder surface, - the negative half-space is unbounded in the x- direction and the - positive half-space is unbounded in all directions. To represent - infinity, numpy.inf is used. - - Parameters - ---------- - side : {'+', '-'} - Indicates the negative or positive half-space - - Returns - ------- - numpy.ndarray - Lower-left coordinates of the axis-aligned bounding box for the - desired half-space - numpy.ndarray - Upper-right coordinates of the axis-aligned bounding box for the - desired half-space - - """ - if side == '-': return (np.array([-np.inf, self.y0 - self.r, self.z0 - self.r]), np.array([np.inf, self.y0 + self.r, self.z0 + self.r])) @@ -1057,57 +1341,21 @@ class XCylinder(Cylinder): np.array([np.inf, np.inf, np.inf])) def evaluate(self, point): - """Evaluate the surface equation at a given point. - - Parameters - ---------- - point : 3-tuple of float - The Cartesian coordinates, :math:`(x',y',z')`, at which the surface - equation should be evaluated. - - Returns - ------- - float - :math:`(y' - y_0)^2 + (z' - z_0)^2 - r^2` - - """ y = point[1] - self.y0 z = point[2] - self.z0 - return y**2 + z**2 - self.r**2 - - def translate(self, vector): - """Translate surface in given direction - - Parameters - ---------- - vector : iterable of float - Direction in which surface should be translated - - Returns - ------- - openmc.XCylinder - Translated surface - - """ - vx, vy, vz = vector - if vy == 0.0 and vz == 0.0: - return self - else: - y0 = self.y0 + vy - z0 = self.z0 + vz - return type(self)(y0=y0, z0=z0, r=self.r) + return y*y + z*z - self.r**2 -class YCylinder(Cylinder): +class YCylinder(QuadricMixin, Surface): """An infinite cylinder whose length is parallel to the y-axis of the form :math:`(x - x_0)^2 + (z - z_0)^2 = r^2`. Parameters ---------- x0 : float, optional - x-coordinate of the center of the cylinder. Defaults to 0. + x-coordinate for the origin of the Cylinder. Defaults to 0 z0 : float, optional - z-coordinate of the center of the cylinder. Defaults to 0. + z-coordinate for the origin of the Cylinder. Defaults to 0 r : float, optional Radius of the cylinder. Defaults to 1. boundary_type : {'transmission, 'vacuum', 'reflective', 'white'}, optional @@ -1124,9 +1372,11 @@ class YCylinder(Cylinder): Attributes ---------- x0 : float - x-coordinate of the center of the cylinder + x-coordinate for the origin of the Cylinder z0 : float - z-coordinate of the center of the cylinder + z-coordinate for the origin of the Cylinder + r : float + Radius of the cylinder boundary_type : {'transmission, 'vacuum', 'reflective', 'white'} Boundary condition that defines the behavior for particles hitting the surface. @@ -1144,14 +1394,16 @@ class YCylinder(Cylinder): _type = 'y-cylinder' _coeff_keys = ('x0', 'z0', 'r') - def __init__(self, x0=0., z0=0., r=1., boundary_type='transmission', - name='', surface_id=None, *, R=None): + def __init__(self, x0=0., z0=0., r=1., **kwargs): + R = kwargs.pop('R', None) if R is not None: - warn(_WARNING_UPPER.format(type(self).__name__, 'r', 'R'), FutureWarning) + warn(_WARNING_UPPER.format(type(self).__name__, 'r', 'R'), + FutureWarning) r = R - super().__init__(r, boundary_type, name, surface_id) - self.x0 = x0 - self.z0 = z0 + super().__init__(**kwargs) + + for key, val in zip(self._coeff_keys, (x0, z0, r)): + setattr(self, key, val) @property def x0(self): @@ -1161,6 +1413,26 @@ class YCylinder(Cylinder): def z0(self): return self.coefficients['z0'] + @property + def r(self): + return self.coefficients['r'] + + @property + def y0(self): + return 0. + + @property + def dx(self): + return 0. + + @property + def dy(self): + return 1. + + @property + def dz(self): + return 0. + @x0.setter def x0(self, x0): check_type('x0 coefficient', x0, Real) @@ -1171,31 +1443,21 @@ class YCylinder(Cylinder): check_type('z0 coefficient', z0, Real) self._coefficients['z0'] = z0 + @r.setter + def r(self, r): + check_type('r coefficient', r, Real) + self._coefficients['r'] = r + + def _get_base_coeffs(self): + x0, z0, r = self.x0, self.z0, self.r + + b = d = e = f = h = 0. + a = c = 1. + g, j, k = -2*x0, -2*z0, x0*x0 + z0*z0 - r*r + + return (a, b, c, d, e, f, g, h, j, k) + def bounding_box(self, side): - """Determine an axis-aligned bounding box. - - An axis-aligned bounding box for surface half-spaces is represented by - its lower-left and upper-right coordinates. For the y-cylinder surface, - the negative half-space is unbounded in the y- direction and the - positive half-space is unbounded in all directions. To represent - infinity, numpy.inf is used. - - Parameters - ---------- - side : {'+', '-'} - Indicates the negative or positive half-space - - Returns - ------- - numpy.ndarray - Lower-left coordinates of the axis-aligned bounding box for the - desired half-space - numpy.ndarray - Upper-right coordinates of the axis-aligned bounding box for the - desired half-space - - """ - if side == '-': return (np.array([self.x0 - self.r, -np.inf, self.z0 - self.r]), np.array([self.x0 + self.r, np.inf, self.z0 + self.r])) @@ -1204,76 +1466,42 @@ class YCylinder(Cylinder): np.array([np.inf, np.inf, np.inf])) def evaluate(self, point): - """Evaluate the surface equation at a given point. - - Parameters - ---------- - point : 3-tuple of float - The Cartesian coordinates, :math:`(x',y',z')`, at which the surface - equation should be evaluated. - - Returns - ------- - float - :math:`(x' - x_0)^2 + (z' - z_0)^2 - r^2` - - """ x = point[0] - self.x0 z = point[2] - self.z0 - return x**2 + z**2 - self.r**2 - - def translate(self, vector): - """Translate surface in given direction - - Parameters - ---------- - vector : iterable of float - Direction in which surface should be translated - - Returns - ------- - openmc.YCylinder - Translated surface - - """ - vx, vy, vz = vector - if vx == 0.0 and vz == 0.0: - return self - else: - x0 = self.x0 + vx - z0 = self.z0 + vz - return type(self)(x0=x0, z0=z0, r=self.r) + return x*x + z*z - self.r**2 -class ZCylinder(Cylinder): +class ZCylinder(QuadricMixin, Surface): """An infinite cylinder whose length is parallel to the z-axis of the form :math:`(x - x_0)^2 + (y - y_0)^2 = r^2`. Parameters ---------- - surface_id : int, optional - Unique identifier for the surface. If not specified, an identifier will - automatically be assigned. + x0 : float, optional + x-coordinate for the origin of the Cylinder. Defaults to 0 + y0 : float, optional + y-coordinate for the origin of the Cylinder. Defaults to 0 + r : float, optional + Radius of the cylinder. Defaults to 1. boundary_type : {'transmission, 'vacuum', 'reflective', 'white'}, optional Boundary condition that defines the behavior for particles hitting the surface. Defaults to transmissive boundary condition where particles freely pass through the surface. - x0 : float, optional - x-coordinate of the center of the cylinder. Defaults to 0. - y0 : float, optional - y-coordinate of the center of the cylinder. Defaults to 0. - r : float, optional - Radius of the cylinder. Defaults to 1. name : str, optional Name of the cylinder. If not specified, the name will be the empty string. + surface_id : int, optional + Unique identifier for the surface. If not specified, an identifier will + automatically be assigned. Attributes ---------- x0 : float - x-coordinate of the center of the cylinder + x-coordinate for the origin of the Cylinder y0 : float - y-coordinate of the center of the cylinder + y-coordinate for the origin of the Cylinder + r : float + Radius of the cylinder boundary_type : {'transmission, 'vacuum', 'reflective', 'white'} Boundary condition that defines the behavior for particles hitting the surface. @@ -1291,14 +1519,16 @@ class ZCylinder(Cylinder): _type = 'z-cylinder' _coeff_keys = ('x0', 'y0', 'r') - def __init__(self, x0=0., y0=0., r=1., boundary_type='transmission', - name='', surface_id=None, *, R=None): + def __init__(self, x0=0., y0=0., r=1., **kwargs): + R = kwargs.pop('R', None) if R is not None: - warn(_WARNING_UPPER.format(type(self).__name__, 'r', 'R'), FutureWarning) + warn(_WARNING_UPPER.format(type(self).__name__, 'r', 'R'), + FutureWarning) r = R - super().__init__(r, boundary_type, name, surface_id) - self.x0 = x0 - self.y0 = y0 + super().__init__(**kwargs) + + for key, val in zip(self._coeff_keys, (x0, y0, r)): + setattr(self, key, val) @property def x0(self): @@ -1308,6 +1538,26 @@ class ZCylinder(Cylinder): def y0(self): return self.coefficients['y0'] + @property + def r(self): + return self.coefficients['r'] + + @property + def z0(self): + return 0. + + @property + def dx(self): + return 0. + + @property + def dy(self): + return 0. + + @property + def dz(self): + return 1. + @x0.setter def x0(self, x0): check_type('x0 coefficient', x0, Real) @@ -1318,31 +1568,21 @@ class ZCylinder(Cylinder): check_type('y0 coefficient', y0, Real) self._coefficients['y0'] = y0 + @r.setter + def r(self, r): + check_type('r coefficient', r, Real) + self._coefficients['r'] = r + + def _get_base_coeffs(self): + x0, y0, r = self.x0, self.y0, self.r + + c = d = e = f = j = 0. + a = b = 1. + g, h, k = -2*x0, -2*y0, x0*x0 + y0*y0 - r*r + + return (a, b, c, d, e, f, g, h, j, k) + def bounding_box(self, side): - """Determine an axis-aligned bounding box. - - An axis-aligned bounding box for surface half-spaces is represented by - its lower-left and upper-right coordinates. For the z-cylinder surface, - the negative half-space is unbounded in the z- direction and the - positive half-space is unbounded in all directions. To represent - infinity, numpy.inf is used. - - Parameters - ---------- - side : {'+', '-'} - Indicates the negative or positive half-space - - Returns - ------- - numpy.ndarray - Lower-left coordinates of the axis-aligned bounding box for the - desired half-space - numpy.ndarray - Upper-right coordinates of the axis-aligned bounding box for the - desired half-space - - """ - if side == '-': return (np.array([self.x0 - self.r, self.y0 - self.r, -np.inf]), np.array([self.x0 + self.r, self.y0 + self.r, np.inf])) @@ -1351,48 +1591,12 @@ class ZCylinder(Cylinder): np.array([np.inf, np.inf, np.inf])) def evaluate(self, point): - """Evaluate the surface equation at a given point. - - Parameters - ---------- - point : 3-tuple of float - The Cartesian coordinates, :math:`(x',y',z')`, at which the surface - equation should be evaluated. - - Returns - ------- - float - :math:`(x' - x_0)^2 + (y' - y_0)^2 - r^2` - - """ x = point[0] - self.x0 y = point[1] - self.y0 - return x**2 + y**2 - self.r**2 - - def translate(self, vector): - """Translate surface in given direction - - Parameters - ---------- - vector : iterable of float - Direction in which surface should be translated - - Returns - ------- - openmc.ZCylinder - Translated surface - - """ - vx, vy, vz = vector - if vx == 0.0 and vy == 0.0: - return self - else: - x0 = self.x0 + vx - y0 = self.y0 + vy - return type(self)(x0=x0, y0=y0, r=self.r) + return x*x + y*y - self.r**2 -class Sphere(Surface): +class Sphere(QuadricMixin, Surface): """A sphere of the form :math:`(x - x_0)^2 + (y - y_0)^2 + (z - z_0)^2 = r^2`. Parameters @@ -1442,16 +1646,16 @@ class Sphere(Surface): _type = 'sphere' _coeff_keys = ('x0', 'y0', 'z0', 'r') - def __init__(self, x0=0., y0=0., z0=0., r=1., boundary_type='transmission', - name='', surface_id=None, *, R=None): + def __init__(self, x0=0., y0=0., z0=0., r=1., **kwargs): + R = kwargs.pop('R', None) if R is not None: - warn(_WARNING_UPPER.format(type(self).__name__, 'r', 'R'), FutureWarning) + warn(_WARNING_UPPER.format(type(self).__name__, 'r', 'R'), + FutureWarning) r = R - super().__init__(surface_id, boundary_type, name=name) - self.x0 = x0 - self.y0 = y0 - self.z0 = z0 - self.r = r + super().__init__(**kwargs) + + for key, val in zip(self._coeff_keys, (x0, y0, z0, r)): + setattr(self, key, val) @property def x0(self): @@ -1489,30 +1693,16 @@ class Sphere(Surface): check_type('r coefficient', r, Real) self._coefficients['r'] = r + def _get_base_coeffs(self): + x0, y0, z0, r = self.x0, self.y0, self.z0, self.r + a = b = c = 1. + d = e = f = 0. + g, h, j = -2*x0, -2*y0, -2*z0 + k = x0*x0 + y0*y0 + z0*z0 - r*r + + return (a, b, c, d, e, f, g, h, j, k) + def bounding_box(self, side): - """Determine an axis-aligned bounding box. - - An axis-aligned bounding box for surface half-spaces is represented by - its lower-left and upper-right coordinates. The positive half-space of a - sphere is unbounded in all directions. To represent infinity, numpy.inf - is used. - - Parameters - ---------- - side : {'+', '-'} - Indicates the negative or positive half-space - - Returns - ------- - numpy.ndarray - Lower-left coordinates of the axis-aligned bounding box for the - desired half-space - numpy.ndarray - Upper-right coordinates of the axis-aligned bounding box for the - desired half-space - - """ - if side == '-': return (np.array([self.x0 - self.r, self.y0 - self.r, self.z0 - self.r]), @@ -1523,50 +1713,13 @@ class Sphere(Surface): np.array([np.inf, np.inf, np.inf])) def evaluate(self, point): - """Evaluate the surface equation at a given point. - - Parameters - ---------- - point : 3-tuple of float - The Cartesian coordinates, :math:`(x',y',z')`, at which the surface - equation should be evaluated. - - Returns - ------- - float - :math:`(x' - x_0)^2 + (y' - y_0)^2 + (z' - z_0)^2 - r^2` - - """ x = point[0] - self.x0 y = point[1] - self.y0 z = point[2] - self.z0 - return x**2 + y**2 + z**2 - self.r**2 - - def translate(self, vector): - """Translate surface in given direction - - Parameters - ---------- - vector : iterable of float - Direction in which surface should be translated - - Returns - ------- - openmc.Sphere - Translated surface - - """ - vx, vy, vz = vector - if vx == 0.0 and vy == 0.0 and vz == 0.0: - return self - else: - x0 = self.x0 + vx - y0 = self.y0 + vy - z0 = self.z0 + vz - return type(self)(x0=x0, y0=y0, z0=z0, r=self.r) + return x*x + y*y + z*z - self.r**2 -class Cone(Surface): +class Cone(QuadricMixin, Surface): """A conical surface parallel to the x-, y-, or z-axis. Parameters @@ -1579,6 +1732,15 @@ class Cone(Surface): z-coordinate of the apex. Defaults to 0. r2 : float, optional Parameter related to the aperature. Defaults to 1. + dx : float, optional + x-component of the vector representing the axis of the cone. + Defaults to 0. + dy : float, optional + y-component of the vector representing the axis of the cone. + Defaults to 0. + dz : float, optional + z-component of the vector representing the axis of the cone. + Defaults to 1. surface_id : int, optional Unique identifier for the surface. If not specified, an identifier will automatically be assigned. @@ -1599,6 +1761,12 @@ class Cone(Surface): z-coordinate of the apex r2 : float Parameter related to the aperature + dx : float + x-component of the vector representing the axis of the cone. + dy : float + y-component of the vector representing the axis of the cone. + dz : float + z-component of the vector representing the axis of the cone. boundary_type : {'transmission, 'vacuum', 'reflective', 'white'} Boundary condition that defines the behavior for particles hitting the surface. @@ -1613,18 +1781,25 @@ class Cone(Surface): """ - _coeff_keys = ('x0', 'y0', 'z0', 'r2') + _type = 'cone' + _coeff_keys = ('x0', 'y0', 'z0', 'r2', 'dx', 'dy', 'dz') - def __init__(self, x0=0., y0=0., z0=0., r2=1., boundary_type='transmission', - name='', surface_id=None, *, R2=None): + def __init__(self, x0=0., y0=0., z0=0., r2=1., dx=0., dy=0., dz=1., **kwargs): + R2 = kwargs.pop('R2', None) if R2 is not None: - warn(_WARNING_UPPER.format(type(self).__name__, 'r2', 'R2'), FutureWarning) + warn(_WARNING_UPPER.format(type(self).__name__, 'r2', 'R2'), + FutureWarning) r2 = R2 - super().__init__(surface_id, boundary_type, name=name) - self.x0 = x0 - self.y0 = y0 - self.z0 = z0 - self.r2 = r2 + super().__init__(**kwargs) + + for key, val in zip(self._coeff_keys, (x0, y0, z0, r2, dx, dy, dz)): + setattr(self, key, val) + + @classmethod + def __subclasshook__(cls, c): + if cls is Cone and c in (XCone, YCone, ZCone): + return True + return NotImplemented @property def x0(self): @@ -1642,6 +1817,18 @@ class Cone(Surface): def r2(self): return self.coefficients['r2'] + @property + def dx(self): + return self.coefficients['dx'] + + @property + def dy(self): + return self.coefficients['dy'] + + @property + def dz(self): + return self.coefficients['dz'] + @x0.setter def x0(self, x0): check_type('x0 coefficient', x0, Real) @@ -1662,31 +1849,72 @@ class Cone(Surface): check_type('r^2 coefficient', r2, Real) self._coefficients['r2'] = r2 - def translate(self, vector): - """Translate surface in given direction + @dx.setter + def dx(self, dx): + check_type('dx coefficient', dx, Real) + self._coefficients['dx'] = dx - Parameters - ---------- - vector : iterable of float - Direction in which surface should be translated + @dy.setter + def dy(self, dy): + check_type('dy coefficient', dy, Real) + self._coefficients['dy'] = dy + + @dz.setter + def dz(self, dz): + check_type('dz coefficient', dz, Real) + self._coefficients['dz'] = dz + + def _get_base_coeffs(self): + # The equation for a general cone with vertex at point p = (x0, y0, z0) + # and axis specified by the unit vector d = (dx, dy, dz) and opening + # half angle theta can be described by the equation + # + # (d*(r - p))^2 - (r - p)*(r - p)cos^2(theta) = 0 + # + # where * is the dot product and the vector r is the evaulation point + # r = (x, y, z) + # + # The argument r2 for cones is actually tan^2(theta) so that + # cos^2(theta) = 1 / (1 + r2) + + x0, y0, z0 = self._origin + dx, dy, dz = self._axis + cos2 = 1 / (1 + self.r2) + + a = cos2 - dx*dx + b = cos2 - dy*dy + c = cos2 - dz*dz + d = -2*dx*dy + e = -2*dy*dz + f = -2*dx*dz + g = 2*(dx*(dy*y0 + dz*z0) - a*x0) + h = 2*(dy*(dx*x0 + dz*z0) - b*y0) + j = 2*(dz*(dx*x0 + dy*y0) - c*z0) + k = a*x0*x0 + b*y0*y0 + c*z0*z0 - 2*(dx*dy*x0*y0 + dy*dz*y0*z0 + + dx*dz*x0*z0) + + return (a, b, c, d, e, f, g, h, j, k) + + def to_xml_element(self): + """Return XML representation of the surface Returns ------- - openmc.Cone - Translated surface + element : xml.etree.ElementTree.Element + XML element containing source data """ - vx, vy, vz = vector - if vx == 0.0 and vy == 0.0 and vz == 0.0: - return self - else: - x0 = self.x0 + vx - y0 = self.y0 + vy - z0 = self.z0 + vz - return type(self)(x0=x0, y0=y0, z0=z0, r2=self.r2) + # This method overrides Surface.to_xml_element to generate a Quadric + # since the C++ layer doesn't support Cones right now + with catch_warnings(): + simplefilter('ignore', IDWarning) + kwargs = {'boundary_type': self.boundary_type, 'name': self.name, + 'surface_id': self.id} + quad_rep = Quadric(*self._get_base_coeffs(), **kwargs) + return quad_rep.to_xml_element() -class XCone(Cone): +class XCone(QuadricMixin, Surface): """A cone parallel to the x-axis of the form :math:`(y - y_0)^2 + (z - z_0)^2 = r^2 (x - x_0)^2`. @@ -1735,29 +1963,86 @@ class XCone(Cone): """ _type = 'x-cone' + _coeff_keys = ('x0', 'y0', 'z0', 'r2') + + def __init__(self, x0=0., y0=0., z0=0., r2=1., **kwargs): + R2 = kwargs.pop('R2', None) + if R2 is not None: + warn(_WARNING_UPPER.format(type(self).__name__, 'r2', 'R2'), + FutureWarning) + r2 = R2 + super().__init__(**kwargs) + + for key, val in zip(self._coeff_keys, (x0, y0, z0, r2)): + setattr(self, key, val) + + @property + def x0(self): + return self.coefficients['x0'] + + @property + def y0(self): + return self.coefficients['y0'] + + @property + def z0(self): + return self.coefficients['z0'] + + @property + def r2(self): + return self.coefficients['r2'] + + @property + def dx(self): + return 1. + + @property + def dy(self): + return 0. + + @property + def dz(self): + return 0. + + @x0.setter + def x0(self, x0): + check_type('x0 coefficient', x0, Real) + self._coefficients['x0'] = x0 + + @y0.setter + def y0(self, y0): + check_type('y0 coefficient', y0, Real) + self._coefficients['y0'] = y0 + + @z0.setter + def z0(self, z0): + check_type('z0 coefficient', z0, Real) + self._coefficients['z0'] = z0 + + @r2.setter + def r2(self, r2): + check_type('r^2 coefficient', r2, Real) + self._coefficients['r2'] = r2 + + def _get_base_coeffs(self): + x0, y0, z0, r2 = self.x0, self.y0, self.z0, self.r2 + + a = -r2 + b = c = 1. + d = e = f = 0. + g, h, j = 2*x0*r2, -2*y0, -2*z0 + k = y0*y0 + z0*z0 - r2*x0*x0 + + return (a, b, c, d, e, f, g, h, j, k) def evaluate(self, point): - """Evaluate the surface equation at a given point. - - Parameters - ---------- - point : 3-tuple of float - The Cartesian coordinates, :math:`(x',y',z')`, at which the surface - equation should be evaluated. - - Returns - ------- - float - :math:`(y' - y_0)^2 + (z' - z_0)^2 - r^2(x' - x_0)^2` - - """ x = point[0] - self.x0 y = point[1] - self.y0 z = point[2] - self.z0 - return y**2 + z**2 - self.r2*x**2 + return y*y + z*z - self.r2*x*x -class YCone(Cone): +class YCone(QuadricMixin, Surface): """A cone parallel to the y-axis of the form :math:`(x - x_0)^2 + (z - z_0)^2 = r^2 (y - y_0)^2`. @@ -1806,29 +2091,86 @@ class YCone(Cone): """ _type = 'y-cone' + _coeff_keys = ('x0', 'y0', 'z0', 'r2') + + def __init__(self, x0=0., y0=0., z0=0., r2=1., **kwargs): + R2 = kwargs.pop('R2', None) + if R2 is not None: + warn(_WARNING_UPPER.format(type(self).__name__, 'r2', 'R2'), + FutureWarning) + r2 = R2 + super().__init__(**kwargs) + + for key, val in zip(self._coeff_keys, (x0, y0, z0, r2)): + setattr(self, key, val) + + @property + def x0(self): + return self.coefficients['x0'] + + @property + def y0(self): + return self.coefficients['y0'] + + @property + def z0(self): + return self.coefficients['z0'] + + @property + def r2(self): + return self.coefficients['r2'] + + @property + def dx(self): + return 0. + + @property + def dy(self): + return 1. + + @property + def dz(self): + return 0. + + @x0.setter + def x0(self, x0): + check_type('x0 coefficient', x0, Real) + self._coefficients['x0'] = x0 + + @y0.setter + def y0(self, y0): + check_type('y0 coefficient', y0, Real) + self._coefficients['y0'] = y0 + + @z0.setter + def z0(self, z0): + check_type('z0 coefficient', z0, Real) + self._coefficients['z0'] = z0 + + @r2.setter + def r2(self, r2): + check_type('r^2 coefficient', r2, Real) + self._coefficients['r2'] = r2 + + def _get_base_coeffs(self): + x0, y0, z0, r2 = self.x0, self.y0, self.z0, self.r2 + + b = -r2 + a = c = 1. + d = e = f = 0. + g, h, j = -2*x0, 2*y0*r2, -2*z0 + k = x0*x0 + z0*z0 - r2*y0*y0 + + return (a, b, c, d, e, f, g, h, j, k) def evaluate(self, point): - """Evaluate the surface equation at a given point. - - Parameters - ---------- - point : 3-tuple of float - The Cartesian coordinates, :math:`(x',y',z')`, at which the surface - equation should be evaluated. - - Returns - ------- - float - :math:`(x' - x_0)^2 + (z' - z_0)^2 - r^2(y' - y_0)^2` - - """ x = point[0] - self.x0 y = point[1] - self.y0 z = point[2] - self.z0 - return x**2 + z**2 - self.r2*y**2 + return x*x + z*z - self.r2*y*y -class ZCone(Cone): +class ZCone(QuadricMixin, Surface): """A cone parallel to the x-axis of the form :math:`(x - x_0)^2 + (y - y_0)^2 = r^2 (z - z_0)^2`. @@ -1877,29 +2219,86 @@ class ZCone(Cone): """ _type = 'z-cone' + _coeff_keys = ('x0', 'y0', 'z0', 'r2') + + def __init__(self, x0=0., y0=0., z0=0., r2=1., **kwargs): + R2 = kwargs.pop('R2', None) + if R2 is not None: + warn(_WARNING_UPPER.format(type(self).__name__, 'r2', 'R2'), + FutureWarning) + r2 = R2 + super().__init__(**kwargs) + + for key, val in zip(self._coeff_keys, (x0, y0, z0, r2)): + setattr(self, key, val) + + @property + def x0(self): + return self.coefficients['x0'] + + @property + def y0(self): + return self.coefficients['y0'] + + @property + def z0(self): + return self.coefficients['z0'] + + @property + def r2(self): + return self.coefficients['r2'] + + @property + def dx(self): + return 0. + + @property + def dy(self): + return 0. + + @property + def dz(self): + return 1. + + @x0.setter + def x0(self, x0): + check_type('x0 coefficient', x0, Real) + self._coefficients['x0'] = x0 + + @y0.setter + def y0(self, y0): + check_type('y0 coefficient', y0, Real) + self._coefficients['y0'] = y0 + + @z0.setter + def z0(self, z0): + check_type('z0 coefficient', z0, Real) + self._coefficients['z0'] = z0 + + @r2.setter + def r2(self, r2): + check_type('r^2 coefficient', r2, Real) + self._coefficients['r2'] = r2 + + def _get_base_coeffs(self): + x0, y0, z0, r2 = self.x0, self.y0, self.z0, self.r2 + + c = -r2 + a = b = 1. + d = e = f = 0. + g, h, j = -2*x0, -2*y0, 2*z0*r2 + k = x0*x0 + y0*y0 - r2*z0*z0 + + return (a, b, c, d, e, f, g, h, j, k) def evaluate(self, point): - """Evaluate the surface equation at a given point. - - Parameters - ---------- - point : 3-tuple of float - The Cartesian coordinates, :math:`(x',y',z')`, at which the surface - equation should be evaluated. - - Returns - ------- - float - :math:`(x' - x_0)^2 + (y' - y_0)^2 - r^2(z' - z_0)^2` - - """ x = point[0] - self.x0 y = point[1] - self.y0 z = point[2] - self.z0 - return x**2 + y**2 - self.r2*z**2 + return x*x + y*y - self.r2*z*z -class Quadric(Surface): +class Quadric(QuadricMixin, Surface): """A surface of the form :math:`Ax^2 + By^2 + Cz^2 + Dxy + Eyz + Fxz + Gx + Hy + Jz + K = 0`. @@ -1939,18 +2338,12 @@ class Quadric(Surface): _coeff_keys = ('a', 'b', 'c', 'd', 'e', 'f', 'g', 'h', 'j', 'k') def __init__(self, a=0., b=0., c=0., d=0., e=0., f=0., g=0., h=0., j=0., - k=0., boundary_type='transmission', name='', surface_id=None): - super().__init__(surface_id, boundary_type, name=name) - self.a = a - self.b = b - self.c = c - self.d = d - self.e = e - self.f = f - self.g = g - self.h = h - self.j = j - self.k = k + k=0., **kwargs): + + super().__init__(**kwargs) + + for key, val in zip(self._coeff_keys, (a, b, c, d, e, f, g, h, j, k)): + setattr(self, key, val) @property def a(self): @@ -2042,50 +2435,8 @@ class Quadric(Surface): check_type('k coefficient', k, Real) self._coefficients['k'] = k - def evaluate(self, point): - """Evaluate the surface equation at a given point. - - Parameters - ---------- - point : 3-tuple of float - The Cartesian coordinates, :math:`(x',y',z')`, at which the surface - equation should be evaluated. - - Returns - ------- - float - :math:`Ax'^2 + By'^2 + Cz'^2 + Dx'y' + Ey'z' + Fx'z' + Gx' + Hy' + - Jz' + K = 0` - - """ - x, y, z = point - return x*(self.a*x + self.d*y + self.g) + \ - y*(self.b*y + self.e*z + self.h) + \ - z*(self.c*z + self.f*x + self.j) + self.k - - def translate(self, vector): - """Translate surface in given direction - - Parameters - ---------- - vector : iterable of float - Direction in which surface should be translated - - Returns - ------- - openmc.Quadric - Translated surface - - """ - vx, vy, vz = vector - a, b, c, d, e, f, g, h, j, k = (getattr(self, key) for key in - self._coeff_keys) - k = (k + vx*vx + vy*vy + vz*vz + d*vx*vy + e*vy*vz + f*vx*vz - - g*vx - h*vy - j*vz) - g = g - 2*a*vx - d*vy - f*vz - h = h - 2*b*vy - d*vx - e*vz - j = j - 2*c*vz - e*vy - f*vx - return type(self)(a=a, b=b, c=c, d=d, e=e, f=f, g=g, h=h, j=j, k=k) + def _get_base_coeffs(self): + return tuple(getattr(self, c) for c in self._coeff_keys) class Halfspace(Region): @@ -2274,3 +2625,5 @@ class Halfspace(Region): # Return translated surface return type(self)(memo[key], self.side) + +_SURFACE_CLASSES = {cls._type: cls for cls in Surface.__subclasses__()} diff --git a/setup.py b/setup.py index 99db83148..c193645c3 100755 --- a/setup.py +++ b/setup.py @@ -56,14 +56,13 @@ kwargs = { 'Topic :: Scientific/Engineering' 'Programming Language :: C++', 'Programming Language :: Python :: 3', - 'Programming Language :: Python :: 3.4', 'Programming Language :: Python :: 3.5', 'Programming Language :: Python :: 3.6', 'Programming Language :: Python :: 3.7', ], # Dependencies - 'python_requires': '>=3.4', + 'python_requires': '>=3.5', 'install_requires': [ 'numpy>=1.9', 'h5py', 'scipy', 'ipython', 'matplotlib', 'pandas', 'lxml', 'uncertainties' diff --git a/src/cell.cpp b/src/cell.cpp index cd6199b75..1d175f05d 100644 --- a/src/cell.cpp +++ b/src/cell.cpp @@ -377,7 +377,7 @@ CSGCell::CSGCell(pugi::xml_node cell_node) throw std::runtime_error{"Invalid surface ID " + std::to_string(abs(r)) + " specified in region for cell " + std::to_string(id_) + "."}; } - r = copysign(it->second + 1, r); + r = (r > 0) ? it->second + 1 : -(it->second + 1); } } @@ -536,8 +536,8 @@ CSGCell::to_hdf5(hid_t cell_group) const region_spec << " |"; } else { // Note the off-by-one indexing - region_spec << " " - << copysign(model::surfaces[abs(token)-1]->id_, token); + auto surf_id = model::surfaces[abs(token)-1]->id_; + region_spec << " " << ((token > 0) ? surf_id : -surf_id); } } write_string(group, "region", region_spec.str(), false); diff --git a/src/output.cpp b/src/output.cpp index d98cfd4a1..18b94be9b 100644 --- a/src/output.cpp +++ b/src/output.cpp @@ -192,7 +192,7 @@ extern "C" void print_particle(Particle* p) // Display miscellaneous info. if (p->surface_ != 0) { const Surface& surf {*model::surfaces[std::abs(p->surface_)-1]}; - fmt::print(" Surface = {}\n", std::copysign(surf.id_, p->surface_)); + fmt::print(" Surface = {}\n", (p->surface_ > 0) ? surf.id_ : -surf.id_); } fmt::print(" Weight = {}\n", p->wgt_); if (settings::run_CE) { diff --git a/src/particle.cpp b/src/particle.cpp index d2d8b2e1b..ea0d03cd0 100644 --- a/src/particle.cpp +++ b/src/particle.cpp @@ -1,7 +1,7 @@ #include "openmc/particle.h" #include // copy, min -#include // log, abs, copysign +#include // log, abs #include @@ -539,7 +539,7 @@ Particle::cross_surface() // TODO: off-by-one surface_ = rotational ? surf_p->i_periodic_ + 1 : - std::copysign(surf_p->i_periodic_ + 1, surface_); + ((surface_ > 0) ? surf_p->i_periodic_ + 1 : -(surf_p->i_periodic_ + 1)); // Figure out what cell particle is in now n_coord_ = 1; diff --git a/src/settings.cpp b/src/settings.cpp index 7808a9d0f..614dd8f5b 100644 --- a/src/settings.cpp +++ b/src/settings.cpp @@ -74,6 +74,7 @@ std::string path_input; std::string path_output; std::string path_particle_restart; std::string path_source; +std::string path_source_library; std::string path_sourcepoint; std::string path_statepoint; diff --git a/src/source.cpp b/src/source.cpp index b15a15d23..96175804c 100644 --- a/src/source.cpp +++ b/src/source.cpp @@ -1,7 +1,15 @@ #include "openmc/source.h" +#if defined (__unix__) || (defined (__APPLE__) && defined (__MACH__)) +#define HAS_DYNAMIC_LINKING +#endif + #include // for move +#ifdef HAS_DYNAMIC_LINKING +#include // for dlopen, dlsym, dlclose, dlerror +#endif + #include #include "xtensor/xadapt.hpp" @@ -71,7 +79,12 @@ SourceDistribution::SourceDistribution(pugi::xml_node node) fatal_error(fmt::format("Source file '{}' does not exist.", settings::path_source)); } - + } else if (check_for_node(node, "library")) { + settings::path_source_library = get_node_value(node, "library", false, true); + if (!file_exists(settings::path_source_library)) { + fatal_error(fmt::format("Source library '{}' does not exist.", + settings::path_source_library)); + } } else { // Spatial distribution for external source @@ -237,7 +250,7 @@ void initialize_source() { write_message("Initializing source particles...", 5); - if (settings::path_source != "") { + if (!settings::path_source.empty()) { // Read the source from a binary file instead of sampling from some // assumed source distribution @@ -261,6 +274,27 @@ void initialize_source() // Close file file_close(file_id); + } else if (!settings::path_source_library.empty()) { + + #ifdef HAS_DYNAMIC_LINKING + // Get the source from a library object + write_message(fmt::format("Sampling from library source {}...", + settings::path_source), 6); + + // Open the library + auto source_library = dlopen(settings::path_source_library.c_str(), RTLD_LAZY); + if (!source_library) { + fatal_error("Couldn't open source library " + settings::path_source_library); + } + + // reset errors + dlerror(); + + fill_source_bank_custom_source(); + #else + fatal_error("Custom source libraries have not yet been implemented for " + "non-POSIX systems"); + #endif } else { // Generation source sites from specified distribution in user input @@ -321,10 +355,50 @@ void free_memory_source() model::external_sources.clear(); } +// fill the source bank from the external source +void fill_source_bank_custom_source() +{ +#ifdef HAS_DYNAMIC_LINKING + // Open the library + auto source_library = dlopen(settings::path_source_library.c_str(), RTLD_LAZY); + if (!source_library) { + fatal_error("Couldn't open source library " + settings::path_source_library); + } + + // reset errors + dlerror(); + + // get the function from the library + using sample_t = Particle::Bank (*)(uint64_t* seed); + auto sample_source = reinterpret_cast(dlsym(source_library, "sample_source")); + + // check for any dlsym errors + auto dlsym_error = dlerror(); + if (dlsym_error) { + dlclose(source_library); + fatal_error(fmt::format("Couldn't open the sample_source symbol: {}", dlsym_error)); + } + + // Generation source sites from specified distribution in the + // library source + for (int64_t i = 0; i < simulation::work_per_rank; ++i) { + // initialize random number seed + int64_t id = (simulation::total_gen + overall_generation()) * + settings::n_particles + simulation::work_index[mpi::rank] + i + 1; + uint64_t seed = init_seed(id, STREAM_SOURCE); + + // sample external source distribution + simulation::source_bank[i] = sample_source(&seed); + } + + // release the library + dlclose(source_library); +#endif +} + void fill_source_bank_fixedsource() { - if (settings::path_source.empty()) { - #pragma omp parallel for + if (settings::path_source.empty() && settings::path_source_library.empty()) { for (int64_t i = 0; i < simulation::work_per_rank; ++i) { // initialize random number seed int64_t id = (simulation::total_gen + overall_generation()) * @@ -334,6 +408,8 @@ void fill_source_bank_fixedsource() // sample external source distribution simulation::source_bank[i] = sample_external_source(&seed); } + } else if (settings::path_source.empty() && !settings::path_source_library.empty()) { + fill_source_bank_custom_source(); } } diff --git a/tests/regression_tests/source_dlopen/__init__.py b/tests/regression_tests/source_dlopen/__init__.py new file mode 100644 index 000000000..e69de29bb diff --git a/tests/regression_tests/source_dlopen/inputs_true.dat b/tests/regression_tests/source_dlopen/inputs_true.dat new file mode 100644 index 000000000..23878ac20 --- /dev/null +++ b/tests/regression_tests/source_dlopen/inputs_true.dat @@ -0,0 +1,23 @@ + + + + + + + + + + + + + + + + + + fixed source + 1000 + 10 + 0 + + diff --git a/tests/regression_tests/source_dlopen/results_true.dat b/tests/regression_tests/source_dlopen/results_true.dat new file mode 100644 index 000000000..e69de29bb diff --git a/tests/regression_tests/source_dlopen/source_sampling.cpp b/tests/regression_tests/source_dlopen/source_sampling.cpp new file mode 100644 index 000000000..eaf8b74d9 --- /dev/null +++ b/tests/regression_tests/source_dlopen/source_sampling.cpp @@ -0,0 +1,23 @@ +#include +#include "openmc/random_lcg.h" +#include "openmc/source.h" +#include "openmc/particle.h" + +// you must have external C linkage here otherwise +// dlopen will not find the file +extern "C" openmc::Particle::Bank sample_source(uint64_t *seed) { + openmc::Particle::Bank particle; + // wgt + particle.particle = openmc::Particle::Type::neutron; + particle.wgt = 1.0; + // position + + particle.r.x = 0.; + particle.r.y = 0.; + particle.r.z = 0.; + // angle + particle.u = {1.0, 0.0, 0.0}; + particle.E = 14.08e6; + particle.delayed_group = 0; + return particle; +} diff --git a/tests/regression_tests/source_dlopen/test.py b/tests/regression_tests/source_dlopen/test.py new file mode 100644 index 000000000..12969cf29 --- /dev/null +++ b/tests/regression_tests/source_dlopen/test.py @@ -0,0 +1,73 @@ +from pathlib import Path +import os +import shutil +import subprocess +import textwrap + +import openmc +import pytest + +from tests.testing_harness import PyAPITestHarness + + +@pytest.fixture +def compile_source(request): + """Compile the external source""" + + # Get build directory and write CMakeLists.txt file + openmc_dir = Path(str(request.config.rootdir)) / 'build' + with open('CMakeLists.txt', 'w') as f: + f.write(textwrap.dedent(""" + cmake_minimum_required(VERSION 3.3 FATAL_ERROR) + project(openmc_sources CXX) + add_library(source SHARED source_sampling.cpp) + find_package(OpenMC REQUIRED HINTS {}) + target_link_libraries(source OpenMC::libopenmc) + """.format(openmc_dir))) + + # Create temporary build directory and change to there + local_builddir = Path('build') + local_builddir.mkdir(exist_ok=True) + os.chdir(str(local_builddir)) + + # Run cmake/make to build the shared libary + subprocess.run(['cmake', os.path.pardir], check=True) + subprocess.run(['make'], check=True) + os.chdir(os.path.pardir) + + yield + + # Remove local build directory when test is complete + shutil.rmtree('build') + + +@pytest.fixture +def model(): + model = openmc.model.Model() + natural_lead = openmc.Material(name="natural_lead") + natural_lead.add_element('Pb', 1.0) + natural_lead.set_density('g/cm3', 11.34) + model.materials.append(natural_lead) + + # geometry + surface_sph1 = openmc.Sphere(r=100, boundary_type='vacuum') + cell_1 = openmc.Cell(fill=natural_lead, region=-surface_sph1) + model.geometry = openmc.Geometry([cell_1]) + + # settings + model.settings.batches = 10 + model.settings.inactive = 0 + model.settings.particles = 1000 + model.settings.run_mode = 'fixed source' + + # custom source from shared library + source = openmc.Source() + source.library = 'build/libsource.so' + model.settings.source = source + + return model + + +def test_dlopen_source(compile_source, model): + harness = PyAPITestHarness('statepoint.10.h5', model) + harness.main() diff --git a/tests/unit_tests/test_data_photon.py b/tests/unit_tests/test_data_photon.py index c767d19e6..f7274e9c1 100644 --- a/tests/unit_tests/test_data_photon.py +++ b/tests/unit_tests/test_data_photon.py @@ -2,6 +2,7 @@ from collections.abc import Mapping, Callable import os +from pathlib import Path import numpy as np import pandas as pd @@ -144,3 +145,9 @@ def test_export_to_hdf5(tmpdir, element): element2.bremsstrahlung['electron_energy']).all() # Export to hdf5 again element2.export_to_hdf5(filename, 'w') + +def test_photodat_only(run_in_tmpdir): + endf_dir = Path(os.environ['OPENMC_ENDF_DATA']) + photoatomic_file = endf_dir / 'photoat' / 'photoat-001_H_000.endf' + data = openmc.data.IncidentPhoton.from_endf(photoatomic_file) + data.export_to_hdf5('tmp.h5', 'w') \ No newline at end of file diff --git a/tests/unit_tests/test_deplete_nuclide.py b/tests/unit_tests/test_deplete_nuclide.py index 7d68b3a3d..b8a258df0 100644 --- a/tests/unit_tests/test_deplete_nuclide.py +++ b/tests/unit_tests/test_deplete_nuclide.py @@ -149,12 +149,18 @@ def test_fission_yield_distribution(): # __getitem__ return yields as a view into yield matrix assert orig_yields.yields.base is yield_dist.yield_matrix - # Fission yield feature uses scaled and incremented + # Scale and increment fission yields mod_yields = orig_yields * 2 assert numpy.array_equal(orig_yields.yields * 2, mod_yields.yields) mod_yields += orig_yields assert numpy.array_equal(orig_yields.yields * 3, mod_yields.yields) + mod_yields = 2.0 * orig_yields + assert numpy.array_equal(orig_yields.yields * 2, mod_yields.yields) + + mod_yields = numpy.float64(2.0) * orig_yields + assert numpy.array_equal(orig_yields.yields * 2, mod_yields.yields) + # Failure modes for adding, multiplying yields similar = numpy.empty_like(orig_yields.yields) with pytest.raises(TypeError): diff --git a/tests/unit_tests/test_source.py b/tests/unit_tests/test_source.py index 1c70e159d..d4d17a3da 100644 --- a/tests/unit_tests/test_source.py +++ b/tests/unit_tests/test_source.py @@ -34,3 +34,11 @@ def test_source_file(): elem = src.to_xml_element() assert 'strength' in elem.attrib assert 'file' in elem.attrib + +def test_source_dlopen(): + library = './libsource.so' + src = openmc.Source(library=library) + assert src.library == library + + elem = src.to_xml_element() + assert 'library' in elem.attrib diff --git a/tests/unit_tests/test_surface.py b/tests/unit_tests/test_surface.py index 703010a05..749ae9621 100644 --- a/tests/unit_tests/test_surface.py +++ b/tests/unit_tests/test_surface.py @@ -2,6 +2,7 @@ from functools import partial from random import uniform, seed import numpy as np +import math import openmc import pytest @@ -139,6 +140,51 @@ def test_zplane(): repr(s) +def test_cylinder(): + x0, y0, z0, r = 2, 3, 4, 2 + dx, dy, dz = 1, -1, 1 + s = openmc.Cylinder(x0=x0, y0=y0, z0=z0, dx=dx, dy=dy, dz=dz, r=r) + assert s.x0 == 2 + assert s.y0 == 3 + assert s.z0 == 4 + assert s.dx == 1 + assert s.dy == -1 + assert s.dz == 1 + assert s.r == 2 + + # Check bounding box + assert_infinite_bb(s) + + # evaluate method + # |(p - p1) ⨯ (p - p2)|^2 / |p2 - p1|^2 - r^2 + p1 = s._origin + p2 = p1 + s._axis + perp = np.array((1, -2, 1))*(1 / s._axis) + divisor = np.linalg.norm(p2 - p1) + pin = p1 + 5*s._axis # point inside cylinder + pout = np.array((4., 0., 2.5)) # point outside the cylinder + pon = p1 + s.r*perp / np.linalg.norm(perp) # point on cylinder + for p, fn in zip((pin, pout, pon), (np.less, np.greater, np.isclose)): + c1 = np.linalg.norm(np.cross(p - p1, p - p2)) / divisor + val = c1*c1 - s.r*s.r + p_eval = s.evaluate(p) + assert fn(p_eval, 0.) + assert p_eval == pytest.approx(val) + + # translate method + st = s.translate((1.0, 1.0, 1.0)) + assert st.x0 == s.x0 + 1 + assert st.y0 == s.y0 + 1 + assert st.z0 == s.z0 + 1 + assert st.dx == s.dx + assert st.dy == s.dy + assert st.dz == s.dz + assert st.r == s.r + + # Make sure repr works + repr(s) + + def test_xcylinder(): y, z, r = 3, 5, 2 s = openmc.XCylinder(y0=y, z0=z, r=r) @@ -284,6 +330,58 @@ def cone_common(apex, r2, cls): repr(s) +def test_cone(): + x0, y0, z0, r2 = 2, 3, 4, 4 + dx, dy, dz = 1, -1, 1 + s = openmc.Cone(x0=x0, y0=y0, z0=z0, dx=dx, dy=dy, dz=dz, r2=r2) + assert s.x0 == 2 + assert s.y0 == 3 + assert s.z0 == 4 + assert s.dx == 1 + assert s.dy == -1 + assert s.dz == 1 + assert s.r2 == 4 + + # Check bounding box + assert_infinite_bb(s) + + # evaluate method + # cos^2(theta) * ((p - p1))**2 - (d @ (p - p1))^2 + # The argument r2 for cones is actually tan^2(theta) so that + # cos^2(theta) = 1 / (1 + r2) + # + # This makes the evaluation equation shown below where p is the evaluation + # point (x, y, z) p1 is the apex (origin) of the cone and r2 is related to + # the aperature of the cone as described above + # (p - p1) @ (p - p1) / (1 + r2) - (d @ (p - p1))^2 + # point inside + p1 = s._origin + d = s._axis + perp = np.array((1, -2, 1))*(1 / d) + perp /= np.linalg.norm(perp) + pin = p1 + 5*d # point inside cone + pout = p1 + 3.2*perp # point outside cone + pon = p1 + 3.2*d + 3.2*math.sqrt(s.r2)*perp # point on cone + for p, fn in zip((pin, pout, pon), (np.less, np.greater, np.isclose)): + val = np.sum((p - p1)**2) / (1 + s.r2) - np.sum((d @ (p - p1))**2) + p_eval = s.evaluate(p) + assert fn(p_eval, 0.) + assert p_eval == pytest.approx(val) + + # translate method + st = s.translate((1.0, 1.0, 1.0)) + assert st.x0 == s.x0 + 1 + assert st.y0 == s.y0 + 1 + assert st.z0 == s.z0 + 1 + assert st.dx == s.dx + assert st.dy == s.dy + assert st.dz == s.dz + assert st.r2 == s.r2 + + # Make sure repr works + repr(s) + + def test_xcone(): apex = (10, 0, 0) r2 = 4 @@ -339,7 +437,7 @@ def test_cylinder_from_points(): p1 = np.array([xi(), xi(), xi()]) p2 = np.array([xi(), xi(), xi()]) r = uniform(1.0, 100.0) - s = openmc.model.cylinder_from_points(p1, p2, r) + s = openmc.Cylinder.from_points(p1, p2, r) # Points p1 and p2 need to be inside cylinder assert p1 in -s @@ -369,24 +467,27 @@ def test_cylinder_from_points_axis(): # (x - 3)^2 + (y - 4)^2 = 2^2 # x^2 + y^2 - 6x - 8y + 21 = 0 - s = openmc.model.cylinder_from_points((3., 4., 0.), (3., 4., 1.), 2.) - assert (s.a, s.b, s.c) == pytest.approx((1., 1., 0.)) - assert (s.d, s.e, s.f) == pytest.approx((0., 0., 0.)) - assert (s.g, s.h, s.j) == pytest.approx((-6., -8., 0.)) - assert s.k == pytest.approx(21.) + s = openmc.Cylinder.from_points((3., 4., 0.), (3., 4., 1.), 2.) + a, b, c, d, e, f, g, h, j, k = s._get_base_coeffs() + assert (a, b, c) == pytest.approx((1., 1., 0.)) + assert (d, e, f) == pytest.approx((0., 0., 0.)) + assert (g, h, j) == pytest.approx((-6., -8., 0.)) + assert k == pytest.approx(21.) # (y + 7)^2 + (z - 1)^2 = 3^2 # y^2 + z^2 + 14y - 2z + 41 = 0 - s = openmc.model.cylinder_from_points((0., -7, 1.), (1., -7., 1.), 3.) - assert (s.a, s.b, s.c) == pytest.approx((0., 1., 1.)) - assert (s.d, s.e, s.f) == pytest.approx((0., 0., 0.)) - assert (s.g, s.h, s.j) == pytest.approx((0., 14., -2.)) - assert s.k == 41. + s = openmc.Cylinder.from_points((0., -7, 1.), (1., -7., 1.), 3.) + a, b, c, d, e, f, g, h, j, k = s._get_base_coeffs() + assert (a, b, c) == pytest.approx((0., 1., 1.)) + assert (d, e, f) == pytest.approx((0., 0., 0.)) + assert (g, h, j) == pytest.approx((0., 14., -2.)) + assert k == 41. # (x - 2)^2 + (z - 5)^2 = 4^2 # x^2 + z^2 - 4x - 10z + 13 = 0 - s = openmc.model.cylinder_from_points((2., 0., 5.), (2., 1., 5.), 4.) - assert (s.a, s.b, s.c) == pytest.approx((1., 0., 1.)) - assert (s.d, s.e, s.f) == pytest.approx((0., 0., 0.)) - assert (s.g, s.h, s.j) == pytest.approx((-4., 0., -10.)) - assert s.k == pytest.approx(13.) + s = openmc.Cylinder.from_points((2., 0., 5.), (2., 1., 5.), 4.) + a, b, c, d, e, f, g, h, j, k = s._get_base_coeffs() + assert (a, b, c) == pytest.approx((1., 0., 1.)) + assert (d, e, f) == pytest.approx((0., 0., 0.)) + assert (g, h, j) == pytest.approx((-4., 0., -10.)) + assert k == pytest.approx(13.)