diff --git a/docs/source/io_formats/settings.rst b/docs/source/io_formats/settings.rst index 2e2648237b..83aa0691f4 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 66440a3c37..ec10fea4fa 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 0000000000..9498176944 --- /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 0000000000..b30884f8ca --- /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 0000000000..606c676df8 --- /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 0000000000..f935ed685e --- /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 0000000000..d68122dd75 --- /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 0000000000..7f6f299261 --- /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 2b4203646f..494cfd1a20 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 a177995ea8..4db17f241a 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/deplete/nuclide.py b/openmc/deplete/nuclide.py index 7e352fbf07..0aae5adb72 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 3ee5052402..9999c383ff 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 88c2f86119..7b709321f1 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 92c828d65f..567968b2c8 100644 --- a/openmc/surface.py +++ b/openmc/surface.py @@ -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__()} diff --git a/setup.py b/setup.py index 99db831484..c193645c35 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/settings.cpp b/src/settings.cpp index 7808a9d0f6..614dd8f5b4 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 b15a15d23a..96175804ce 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 0000000000..e69de29bb2 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 0000000000..23878ac207 --- /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 0000000000..e69de29bb2 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 0000000000..eaf8b74d92 --- /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 0000000000..12969cf29d --- /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_deplete_nuclide.py b/tests/unit_tests/test_deplete_nuclide.py index 7d68b3a3d8..b8a258df00 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 1c70e159d2..d4d17a3dab 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 703010a059..749ae96212 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.)