Merge branch 'develop' into cam_WP7_task3

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 <project>_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})

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()

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")

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 <path to openmc>)
+   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

docs/source/usersguide/install.rst

@@ -404,7 +404,7 @@ to install the Python package in :ref:`"editable" mode <devguide_editable>`.
 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

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)

examples/xml/custom_source/geometry.xml

@@ -0,0 +1,15 @@
+<?xml version="1.0"?>
+<geometry>
+
+  <!-- Definition of Cells -->
+  <cell id="1" universe="0" fill="37" region="-2" />
+  <cell id="100" universe="37" material="40" region="-1" />
+  <cell id="101" universe="37" material="41" region="1" />
+  <cell id="2" universe="0" material="41" region="2 -3" />
+
+  <!-- Defition of Surfaces -->
+  <surface id="1" type="z-cylinder" coeffs="0 0 7" />
+  <surface id="2" type="z-cylinder" coeffs="0 0 9" />
+  <surface id="3" type="z-cylinder" coeffs="0 0 11" boundary="vacuum" />
+
+</geometry>

examples/xml/custom_source/materials.xml

@@ -0,0 +1,16 @@
+<?xml version="1.0"?>
+<materials>
+
+  <material id="40">
+    <density value="4.5" units="g/cc" />
+    <nuclide name="U235" ao="1.0" />
+  </material>
+
+  <material id="41">
+    <density value="1.0" units="g/cc" />
+    <nuclide name="H1" ao="2.0" />
+    <nuclide name="O16" ao="1.0" />
+    <sab name="c_H_in_H2O"/>
+  </material>
+
+</materials>

examples/xml/custom_source/settings.xml

@@ -0,0 +1,14 @@
+<?xml version="1.0"?>
+<settings>
+
+  <run_mode>fixed source</run_mode>
+  <batches>10</batches>
+  <inactive>0</inactive>
+  <particles>100000</particles>
+
+  <!-- Starting source -->
+  <source>
+    <library>build/libsource.so</library>
+  </source>
+
+</settings>

examples/xml/custom_source/source_ring.cpp

@@ -0,0 +1,26 @@
+#include <cmath> // 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;
+}

examples/xml/custom_source/tallies.xml

@@ -0,0 +1,17 @@
+<?xml version="1.0"?>
+<tallies>
+
+  <filter id="1" type="cell">
+    <bins>100</bins>
+  </filter>
+
+  <filter id="2" type="energy">
+    <bins>0 20.0e6</bins>
+  </filter>
+
+  <tally id="3">
+    <filters>1 2 </filters>
+    <scores>flux</scores>
+  </tally>
+
+</tallies>

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
 

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

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
 

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

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):

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)

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'

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);

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) {

src/particle.cpp

@@ -1,7 +1,7 @@
 #include "openmc/particle.h"
 
 #include <algorithm> // copy, min
-#include <cmath>     // log, abs, copysign
+#include <cmath>     // log, abs
 
 #include <fmt/core.h>
 
@@ -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;

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;
 

src/source.cpp

@@ -1,7 +1,15 @@
 #include "openmc/source.h"
 
+#if defined (__unix__) || (defined (__APPLE__) && defined (__MACH__))
+#define HAS_DYNAMIC_LINKING
+#endif
+
 #include <algorithm> // for move
 
+#ifdef HAS_DYNAMIC_LINKING
+#include <dlfcn.h> // for dlopen, dlsym, dlclose, dlerror
+#endif
+
 #include <fmt/core.h>
 #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<sample_t>(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();
   }
 }
 

tests/regression_tests/source_dlopen/inputs_true.dat

@@ -0,0 +1,23 @@
+<?xml version='1.0' encoding='utf-8'?>
+<geometry>
+  <cell id="1" material="1" region="-1" universe="1" />
+  <surface boundary="vacuum" coeffs="0.0 0.0 0.0 100" id="1" type="sphere" />
+</geometry>
+<?xml version='1.0' encoding='utf-8'?>
+<materials>
+  <material id="1" name="natural_lead">
+    <density units="g/cm3" value="11.34" />
+    <nuclide ao="0.014" name="Pb204" />
+    <nuclide ao="0.241" name="Pb206" />
+    <nuclide ao="0.221" name="Pb207" />
+    <nuclide ao="0.524" name="Pb208" />
+  </material>
+</materials>
+<?xml version='1.0' encoding='utf-8'?>
+<settings>
+  <run_mode>fixed source</run_mode>
+  <particles>1000</particles>
+  <batches>10</batches>
+  <inactive>0</inactive>
+  <source library="build/libsource.so" strength="1.0" />
+</settings>

tests/regression_tests/source_dlopen/source_sampling.cpp

@@ -0,0 +1,23 @@
+#include <iostream>
+#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;    
+}

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()

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')

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):

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

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.)