crp/OpenMC
1
0
Fork 0
mirror of https://github.com/openmc-dev/openmc.git synced 2026-07-28 14:15:42 -04:00
Code Issues Projects Releases Packages Wiki Activity Actions 5

Merge remote-tracking branch 'upstream/develop' into develop

Browse source
This commit is contained in:
Simon Richards 2020-02-25 20:08:57 +00:00
commit 663b31af43
25 changed files with 821 additions and 183 deletions
Download patch file Download diff file Expand all files Collapse all files

72
docs/source/io_formats/settings.rst
View file

@ -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

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

@ -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

8
examples/xml/custom_source/CMakeLists.txt Normal file
View file

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

15
examples/xml/custom_source/geometry.xml Normal file
View file

@ -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>

16
examples/xml/custom_source/materials.xml Normal file
View file

@ -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>

14
examples/xml/custom_source/settings.xml Normal file
View file

@ -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>

26
examples/xml/custom_source/source_ring.cpp Normal file
View file

@ -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;
}

17
examples/xml/custom_source/tallies.xml Normal file
View file

@ -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>

1
include/openmc/settings.h
View file

@ -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

3
include/openmc/source.h
View file

@ -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

8
openmc/deplete/nuclide.py
View file

@ -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

47
openmc/model/funcs.py
View file

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

38
openmc/source.py
View file

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

381
openmc/surface.py
View file

@ -3,13 +3,14 @@ from collections import OrderedDict
from copy import deepcopy
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']
@ -361,6 +362,59 @@ class PlaneMixin(metaclass=ABCMeta):
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.
@ -509,6 +563,12 @@ class Plane(PlaneMixin, Surface):
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):
return self.coefficients['a']
@ -545,9 +605,6 @@ class Plane(PlaneMixin, Surface):
check_type('D coefficient', d, Real)
self._coefficients['d'] = d
def _get_base_coeffs(self):
return (self.a, self.b, self.c, self.d)
@classmethod
def from_points(cls, p1, p2, p3, **kwargs):
"""Return a plane given three points that pass through it.
@ -638,29 +695,31 @@ class XPlane(PlaneMixin, Surface):
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 _get_base_coeffs(self):
return (1., 0., 0., self.x0)
def bounding_box(self, side):
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):
return point[0] - self.x0
Plane.register(XPlane)
class YPlane(PlaneMixin, Surface):
"""A plane perpendicular to the y axis of the form :math:`y - y_0 = 0`
@ -718,29 +777,31 @@ class YPlane(PlaneMixin, Surface):
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 _get_base_coeffs(self):
return (0., 1., 0., self.y0)
def bounding_box(self, side):
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):
return point[1] - self.y0
Plane.register(YPlane)
class ZPlane(PlaneMixin, Surface):
"""A plane perpendicular to the z axis of the form :math:`z - z_0 = 0`
@ -798,32 +859,43 @@ class ZPlane(PlaneMixin, Surface):
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 _get_base_coeffs(self):
return (0., 0., 1., self.z0)
def bounding_box(self, side):
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]))
def evaluate(self, point):
return point[2] - self.z0
Plane.register(ZPlane)
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.
@ -884,7 +956,7 @@ class QuadricMixin(metaclass=ABCMeta):
"""
x = np.asarray(point)
A, b, c = self.get_Abc()
return np.matmul(x.T, np.matmul(A, x)) + np.matmul(b.T, x) + c
return x.T @ A @ x + b.T @ x + c
def translate(self, vector, inplace=False):
"""Translate surface in given direction
@ -907,17 +979,20 @@ class QuadricMixin(metaclass=ABCMeta):
surf = self if inplace else self.clone()
if set(('x0', 'y0', 'z0')).intersection(set(surf._coeff_keys)):
if hasattr(self, 'x0'):
for vi, xi in zip(vector, ('x0', 'y0', 'z0')):
val = getattr(surf, xi, None)
if val is not None:
val = getattr(surf, xi)
try:
setattr(surf, xi, val + vi)
except AttributeError:
# That attribute is read only i.e x0 for XCylinder
pass
else:
A, bvec, cnst = self.get_Abc()
g, h, j = bvec - 2*np.matmul(vector.T, A)
k = cnst + np.matmul(vector.T, np.matmul(A, vector)) \
- np.matmul(bvec.T, vector)
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)
@ -993,16 +1068,17 @@ class Cylinder(QuadricMixin, Surface):
_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):
raise NotImplementedError('There is no C++ implementation for general '
'Cylinders yet, please use '
'openmc.model.funcs.cylinder_from_points to '
'return a Quadric instance instead for now')
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']
@ -1068,8 +1144,8 @@ class Cylinder(QuadricMixin, Surface):
def _get_base_coeffs(self):
# Get x, y, z coordinates of two points
x1, y1, z1 = self.x0, self.y0, self.z0
x2, y2, z2 = x1 + self.dx, y1 + self.dy, z1 + self.dz
x1, y1, z1 = self._origin
x2, y2, z2 = self._origin + self._axis
r = self.r
# Define intermediate terms
@ -1118,10 +1194,6 @@ class Cylinder(QuadricMixin, Surface):
radius r.
"""
raise NotImplementedError('There is no C++ implementation for general '
'Cylinders yet, please use '
'openmc.model.funcs.cylinder_from_points to '
'return a Quadric instance instead for now')
# Convert to numpy arrays
p1 = np.asarray(p1)
p2 = np.asarray(p2)
@ -1130,6 +1202,24 @@ class Cylinder(QuadricMixin, Surface):
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
@ -1202,6 +1292,22 @@ class XCylinder(QuadricMixin, Surface):
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)
@ -1240,9 +1346,6 @@ class XCylinder(QuadricMixin, Surface):
return y*y + z*z - self.r**2
Cylinder.register(XCylinder)
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`.
@ -1314,6 +1417,22 @@ class YCylinder(QuadricMixin, Surface):
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)
@ -1352,9 +1471,6 @@ class YCylinder(QuadricMixin, Surface):
return x*x + z*z - self.r**2
Cylinder.register(YCylinder)
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`.
@ -1426,6 +1542,22 @@ class ZCylinder(QuadricMixin, Surface):
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)
@ -1464,9 +1596,6 @@ class ZCylinder(QuadricMixin, Surface):
return x*x + y*y - self.r**2
Cylinder.register(ZCylinder)
class Sphere(QuadricMixin, Surface):
"""A sphere of the form :math:`(x - x_0)^2 + (y - y_0)^2 + (z - z_0)^2 = r^2`.
@ -1656,9 +1785,6 @@ class Cone(QuadricMixin, Surface):
_coeff_keys = ('x0', 'y0', 'z0', 'r2', 'dx', 'dy', 'dz')
def __init__(self, x0=0., y0=0., z0=0., r2=1., dx=0., dy=0., dz=1., **kwargs):
raise NotImplementedError('There is no C++ implementation for general '
'Cones yet, this functionality should be '
'added soon.')
R2 = kwargs.pop('R2', None)
if R2 is not None:
warn(_WARNING_UPPER.format(type(self).__name__, 'r2', 'R2'),
@ -1669,6 +1795,12 @@ class Cone(QuadricMixin, Surface):
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):
return self.coefficients['x0']
@ -1745,27 +1877,42 @@ class Cone(QuadricMixin, Surface):
# The argument r2 for cones is actually tan^2(theta) so that
# cos^2(theta) = 1 / (1 + r2)
x0, y0, z0, r2 = self.x0, self.y0, self.z0, self.r2
dx, dy, dz = self.dx, self.dy, self.dz
dnorm = dx*dx + dy*dy + dz*dz
dx /= dnorm
dy /= dnorm
dz /= dnorm
cos2 = 1 / (1 + r2)
x0, y0, z0 = self._origin
dx, dy, dz = self._axis
cos2 = 1 / (1 + self.r2)
a = dx*dx - cos2
b = dy*dy - cos2
c = dz*dz - cos2
d = 2*dx*dy
e = 2*dy*dz
f = 2*dx*dz
g = -2*(dx*dx*x0 + dx*dy*y0 + dx*dz*z0 - cos2)
h = -2*(dy*dy*y0 + dx*dy*x0 + dy*dz*z0 - cos2)
j = -2*(dz*dz*y0 + dx*dz*x0 + dy*dz*y0 - cos2)
k = (dx*x0 + dy*y0 + dz*z0)**2 - cos2*(x0*x0 + y0*y0 + z0*z0)
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
-------
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 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(QuadricMixin, Surface):
"""A cone parallel to the x-axis of the form :math:`(y - y_0)^2 + (z - z_0)^2 =
@ -1845,6 +1992,18 @@ class XCone(QuadricMixin, Surface):
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)
@ -1883,9 +2042,6 @@ class XCone(QuadricMixin, Surface):
return y*y + z*z - self.r2*x*x
Cone.register(XCone)
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`.
@ -1964,6 +2120,18 @@ class YCone(QuadricMixin, Surface):
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)
@ -2002,9 +2170,6 @@ class YCone(QuadricMixin, Surface):
return x*x + z*z - self.r2*y*y
Cone.register(YCone)
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`.
@ -2083,6 +2248,18 @@ class ZCone(QuadricMixin, Surface):
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)
@ -2121,9 +2298,6 @@ class ZCone(QuadricMixin, Surface):
return x*x + y*y - self.r2*z*z
Cone.register(ZCone)
class Quadric(QuadricMixin, Surface):
"""A surface of the form :math:`Ax^2 + By^2 + Cz^2 + Dxy + Eyz + Fxz + Gx + Hy +
Jz + K = 0`.
@ -2452,5 +2626,4 @@ class Halfspace(Region):
# Return translated surface
return type(self)(memo[key], self.side)
_SURFACE_CLASSES = {cls._type: cls for cls in Surface.__subclasses__()}

3
setup.py
View file

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

1
src/settings.cpp
View file

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

84
src/source.cpp
View file

@ -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();
}
}

0
tests/regression_tests/source_dlopen/__init__.py Normal file
View file

23
tests/regression_tests/source_dlopen/inputs_true.dat Normal file
View file

@ -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>

0
tests/regression_tests/source_dlopen/results_true.dat Normal file
View file

23
tests/regression_tests/source_dlopen/source_sampling.cpp Normal file
View file

@ -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;
}

73
tests/regression_tests/source_dlopen/test.py Normal file
View file

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

8
tests/unit_tests/test_deplete_nuclide.py
View file

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

8
tests/unit_tests/test_source.py
View file

@ -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

133
tests/unit_tests/test_surface.py
View file

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