diff --git a/CMakeLists.txt b/CMakeLists.txt index d0c2569c34..fb5da6f2a6 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -364,6 +364,10 @@ list(APPEND libopenmc_SOURCES src/xml_interface.cpp src/xsdata.cpp) +# Add bundled external dependencies +list(APPEND libopenmc_SOURCES + src/external/quartic_solver.cpp) + # For Visual Studio compilers if(MSVC) # Use static library (otherwise explicit symbol portings are needed) diff --git a/CONTRIBUTING.md b/CONTRIBUTING.md index be506f919a..378ca346a8 100644 --- a/CONTRIBUTING.md +++ b/CONTRIBUTING.md @@ -13,9 +13,8 @@ openmc@anl.gov. ## Resources - [GitHub Repository](https://github.com/openmc-dev/openmc) -- [Documentation](http://openmc.readthedocs.io/en/latest) -- [User's Mailing List](openmc-users@googlegroups.com) -- [Developer's Mailing List](openmc-dev@googlegroups.com) +- [Documentation](http://docs.openmc.org/en/latest) +- [Discussion Forum](https://openmc.discourse.group) - [Slack Community](https://openmc.slack.com/signup) (If you don't see your domain listed, contact openmc@anl.gov) @@ -37,10 +36,10 @@ development team will be happy to discuss it. All changes to OpenMC happen through pull requests. For a full overview of the process, see the developer's guide section on [Contributing to -OpenMC](http://openmc.readthedocs.io/en/latest/devguide/contributing.html). +OpenMC](https://docs.openmc.org/en/latest/devguide/contributing.html). ## Code Style Before you run off to make changes to the code, please have a look at our [style -guide](http://openmc.readthedocs.io/en/latest/devguide/styleguide.html), which +guide](https://docs.openmc.org/en/latest/devguide/styleguide.html), which is used when reviewing new contributions. diff --git a/README.md b/README.md index ebd7b050ae..186c1e5e03 100644 --- a/README.md +++ b/README.md @@ -1,6 +1,6 @@ # OpenMC Monte Carlo Particle Transport Code -[![License](https://img.shields.io/github/license/openmc-dev/openmc.svg)](http://openmc.readthedocs.io/en/latest/license.html) +[![License](https://img.shields.io/badge/license-MIT-green)](https://docs.openmc.org/en/latest/license.html) [![GitHub Actions build status (Linux)](https://github.com/openmc-dev/openmc/workflows/CI/badge.svg?branch=develop)](https://github.com/openmc-dev/openmc/actions?query=workflow%3ACI) [![Code Coverage](https://coveralls.io/repos/github/openmc-dev/openmc/badge.svg?branch=develop)](https://coveralls.io/github/openmc-dev/openmc?branch=develop) [![dockerhub-publish-develop-dagmc](https://github.com/openmc-dev/openmc/workflows/dockerhub-publish-develop-dagmc/badge.svg)](https://github.com/openmc-dev/openmc/actions?query=workflow%3Adockerhub-publish-develop-dagmc) @@ -12,14 +12,14 @@ continuous-energy transport code that uses HDF5 format cross sections. The project started under the Computational Reactor Physics Group at MIT. Complete documentation on the usage of OpenMC is hosted on Read the Docs (both -for the [latest release](http://openmc.readthedocs.io/en/stable/) and -[developmental](http://openmc.readthedocs.io/en/latest/) version). If you are +for the [latest release](https://docs.openmc.org/en/stable/) and +[developmental](https://docs.openmc.org/en/latest/) version). If you are interested in the project, or would like to help and contribute, please get in touch on the OpenMC [discussion forum](https://openmc.discourse.group/). ## Installation Detailed [installation -instructions](http://openmc.readthedocs.io/en/stable/usersguide/install.html) +instructions](https://docs.openmc.org/en/stable/usersguide/install.html) can be found in the User's Guide. ## Citing @@ -35,7 +35,7 @@ citing the following publication: ## Troubleshooting If you run into problems compiling, installing, or running OpenMC, first check -the [Troubleshooting section](http://openmc.readthedocs.io/en/stable/usersguide/troubleshoot.html) in +the [Troubleshooting section](https://docs.openmc.org/en/stable/usersguide/troubleshoot.html) in the User's Guide. If you are not able to find a solution to your problem there, please post to the [discussion forum](https://openmc.discourse.group/). @@ -53,4 +53,4 @@ create an Issue on github. ## License OpenMC is distributed under the MIT/X -[license](http://openmc.readthedocs.io/en/stable/license.html). +[license](https://docs.openmc.org/en/stable/license.html). diff --git a/docs/source/io_formats/geometry.rst b/docs/source/io_formats/geometry.rst index b04edca2b7..ff2f6eb74e 100644 --- a/docs/source/io_formats/geometry.rst +++ b/docs/source/io_formats/geometry.rst @@ -26,7 +26,7 @@ Each ```` element can have the following attributes or sub-elements: :type: The type of the surfaces. This can be "x-plane", "y-plane", "z-plane", "plane", "x-cylinder", "y-cylinder", "z-cylinder", "sphere", "x-cone", - "y-cone", "z-cone", or "quadric". + "y-cone", "z-cone", "quadric", "x-torus", "y-torus", or "z-torus". *Default*: None diff --git a/docs/source/io_formats/summary.rst b/docs/source/io_formats/summary.rst index 9b92b30187..c140d98ae3 100644 --- a/docs/source/io_formats/summary.rst +++ b/docs/source/io_formats/summary.rst @@ -55,7 +55,8 @@ The current version of the summary file format is 6.0. :Datasets: - **name** (*char[]*) -- Name of the surface. - **type** (*char[]*) -- Type of the surface. Can be 'x-plane', 'y-plane', 'z-plane', 'plane', 'x-cylinder', 'y-cylinder', - 'z-cylinder', 'sphere', 'x-cone', 'y-cone', 'z-cone', or 'quadric'. + 'z-cylinder', 'sphere', 'x-cone', 'y-cone', 'z-cone', 'quadric', + 'x-torus', 'y-torus', or 'z-torus'. - **coefficients** (*double[]*) -- Array of coefficients that define the surface. See :ref:`surface_element` for what coefficients are defined for each surface type. diff --git a/docs/source/methods/geometry.rst b/docs/source/methods/geometry.rst index 2c7b7ecfb2..2da0c02112 100644 --- a/docs/source/methods/geometry.rst +++ b/docs/source/methods/geometry.rst @@ -118,7 +118,19 @@ to fully define the surface. +----------------------+------------+------------------------------+-------------------------+ | General quadric | quadric | :math:`Ax^2 + By^2 + Cz^2 + | :math:`A \; B \; C \; D | | surface | | Dxy + Eyz + Fxz + Gx + Hy + | \; E \; F \; G \; H \; | - | | | Jz + K` | J \; K` | + | | | Jz + K = 0` | J \; K` | + +----------------------+------------+------------------------------+-------------------------+ + | Torus parallel to the| x-torus | :math:`(x-x_0)^2/B^2+\frac{( | :math:`x_0 \; y_0 \; | + | :math:`x`-axis | | \sqrt{(y-y_0)^2+(z-z_0)^2} - | z_0 \; A \; B \; C` | + | | | A)^2}{C^2} - 1 = 0` | | + +----------------------+------------+------------------------------+-------------------------+ + | Torus parallel to the| y-torus | :math:`(y-y_0)^2/B^2+\frac{( | :math:`x_0 \; y_0 \; | + | :math:`y`-axis | | \sqrt{(x-x_0)^2+(z-z_0)^2} - | z_0 \; A \; B \; C` | + | | | A)^2}{C^2} - 1 = 0` | | + +----------------------+------------+------------------------------+-------------------------+ + | Torus parallel to the| z-torus | :math:`(z-z_0)^2/B^2+\frac{( | :math:`x_0 \; y_0 \; | + | :math:`z`-axis | | \sqrt{(x-x_0)^2+(y-y_0)^2} - | z_0 \; A \; B \; C` | + | | | A)^2}{C^2} - 1 = 0` | | +----------------------+------------+------------------------------+-------------------------+ .. _universes: @@ -437,6 +449,72 @@ Defining the terms we then have the simple quadratic equation :math:`ad^2 + 2kd + c = 0` which can be solved as described in :ref:`cylinder_distance`. +Torus Parallel to an Axis +------------------------- + +The equation for a torus parallel to, for example, the x-axis is + +.. math:: + :label: dist-xtorus-sqrt + + \frac{(x-x_0)^2}{B^2} + \frac{(\sqrt{(y-y_0)^2 + (z-z_0)^2} - A)^2}{C^2} - + 1 = 0. + +First, it needs to be cast into a polynomial form. Rearranging terms, + +.. math:: + :label: dist-xtorus-1 + + (D\bar{x}^2 + \bar{y}^2 + \bar{z}^2 + A^2 - C^2)^2 = 4A^2(\bar{y}^2 + + \bar{z}^2) + +where :math:`D = (C/B)^2`, :math:`\bar{x} = x - x_0`, :math:`\bar{y} = y - y_0`, +and :math:`\bar{z} = z - z_0`. To find the distance to the surface, we thus need +to solve + +.. math:: + :label: dist-xtorus-2 + + (D(\bar{x} + du)^2 + (\bar{y} + dv)^2 + (\bar{z} + dw)^2 + A^2 - C^2)^2 = + 4A^2((\bar{y} + dv)^2 + (\bar{z} + dw)^2). + +Expanding and collecting like powers of :math:`d` yields + +.. math:: + :label: dist-xtorus-3 + + (c_2d^2 + c_1d + c_0)^2 = c_2'd^2 + c_1'd + c_0' + +where + +.. math:: + :label: dist-xtorus-4 + + \begin{align*} + c_2 &= Du^2 + v^2 + w^2 \\ + c_1 &= 2(Du\bar{x} + v\bar{y} + w\bar{z}) \\ + c_0 &= D\bar{x}^2 + \bar{y}^2 + \bar{z}^2 + A^2 - C^2 \\ + c_2' &= 4A^2 (v^2 + w^2) \\ + c_1' &= 8A^2 (v\bar{y} + w\bar{z}) \\ + c_0' &= 4A^2(\bar{y}^2 + \bar{z}^2). + \end{align*} + +Expanding the left-hand side and collecting like powers of :math:`d` on one +side, we obtain + +.. math:: + :label: dist-xtorus-5 + + (c_2^2)d^4 + (2c_1c_2)d^3 + (c_1^2 + 2c_0c_2 - c_2')d^2 + (2c_0c_1 - c_1')d + + (c_0^2 - c_0') = 0. + +The above equation is a fourth-order (quartic) polynomial equation. Although +there is an analytical solution to the general quartic equation, it can be +subject to roundoff errors when evaluated numerically. OpenMC uses an external +`quartic equation solver `_ developed by +Orellana and De Michele that is based on the decomposition of the quartic +polynomial into two quadratics. + .. _find-cell: ---------------------------- @@ -899,6 +977,27 @@ Dxy + Eyz + Fxz + Gx + Hy + Jz + K = 0`. Thus, the gradient to the surface is \nabla f = \left ( \begin{array}{c} 2Ax + Dy + Fz + G \\ 2By + Dx + Ez + H \\ 2Cz + Ey + Fx + J \end{array} \right ). +Torus Parallel to an Axis +------------------------- + +A torus parallel to, for example, the x-axis has the form + +.. math:: + :label: reflection-torus-1 + + f(x,y,z) = \frac{(x-x_0)^2}{B^2} + \frac{(\sqrt{(y-y_0)^2 + (z-z_0)^2} - + A)^2}{C^2} - 1. + +The gradient to the surface is therefore + +.. math:: + :label: reflection-torus-grad + + \nabla f = \left ( \begin{array}{c} 2\bar{x}/B^2 \\ 2\bar{y}(g - A)/(C^2g) + \\ 2\bar{z}(g - A)/(C^2g) \end{array} \right ) + +where :math:`g = \sqrt{\bar{y}^2 + \bar{z}^2}` and, as always, :math:`\bar{x} = +x - x_0`, :math:`\bar{y} = y - y_0`, and :math:`\bar{z} = z - z_0`. .. _white: diff --git a/docs/source/pythonapi/base.rst b/docs/source/pythonapi/base.rst index 342b4c1a64..496201ad77 100644 --- a/docs/source/pythonapi/base.rst +++ b/docs/source/pythonapi/base.rst @@ -80,6 +80,9 @@ Building geometry openmc.YCone openmc.ZCone openmc.Quadric + openmc.XTorus + openmc.YTorus + openmc.ZTorus openmc.Halfspace openmc.Intersection openmc.Union diff --git a/docs/source/usersguide/geometry.rst b/docs/source/usersguide/geometry.rst index 2b6a6fef2f..68ecfce95c 100644 --- a/docs/source/usersguide/geometry.rst +++ b/docs/source/usersguide/geometry.rst @@ -78,10 +78,23 @@ classes are listed in the following table. | Cone parallel to the | :math:`(x-x_0)^2 + (y-y_0)^2 | :class:`openmc.ZCone` | | :math:`z`-axis | - R^2(z-z_0)^2 = 0` | | +----------------------+------------------------------+---------------------------+ - | General quadric | :math:`Ax^2 + By^2 + Cz^2 + | :class:`openmc.Quadric` | + | General quadric | :math:`Ax^2 + By^2 + Cz^2 + | :class:`openmc.Quadric` | | surface | Dxy + Eyz + Fxz + Gx + Hy + | | | | Jz + K = 0` | | +----------------------+------------------------------+---------------------------+ + | Torus parallel to the| :math:`(x-x_0)^2/B^2+\frac{( | :class:`openmc.XTorus` | + | :math:`x`-axis | \sqrt{(y-y_0)^2+(z-z_0)^2} - | | + | | A)^2}{C^2} - 1 = 0` | | + +----------------------+------------------------------+---------------------------+ + | Torus parallel to the| :math:`(y-y_0)^2/B^2+\frac{( | :class:`openmc.YTorus` | + | :math:`y`-axis | \sqrt{(x-x_0)^2+(z-z_0)^2} - | | + | | A)^2}{C^2} - 1 = 0` | | + +----------------------+------------------------------+---------------------------+ + | Torus parallel to the| :math:`(z-z_0)^2/B^2+\frac{( | :class:`openmc.ZTorus` | + | :math:`z`-axis | \sqrt{(x-x_0)^2+(y-y_0)^2} - | | + | | A)^2}{C^2} - 1 = 0` | | + +----------------------+------------------------------+---------------------------+ + Each surface is characterized by several parameters. As one example, the parameters for a sphere are the :math:`x,y,z` coordinates of the center of the diff --git a/include/openmc/external/quartic_solver.h b/include/openmc/external/quartic_solver.h new file mode 100644 index 0000000000..376a014af6 --- /dev/null +++ b/include/openmc/external/quartic_solver.h @@ -0,0 +1,10 @@ +#ifndef OPENMC_EXTERNAL_QUARTIC_SOLVER_H +#define OPENMC_EXTERNAL_QUARTIC_SOLVER_H + +#include + +namespace oqs { +void quartic_solver(double coeff[5], std::complex roots[4]); +} // end namespace oqs + +#endif // OPENMC_EXTERNAL_QUARTIC_SOLVER_H diff --git a/include/openmc/surface.h b/include/openmc/surface.h index 325c313aa4..84a2ed1134 100644 --- a/include/openmc/surface.h +++ b/include/openmc/surface.h @@ -371,6 +371,57 @@ public: double A_, B_, C_, D_, E_, F_, G_, H_, J_, K_; }; +//============================================================================== +//! A toroidal surface described by the quartic torus lies in the x direction +// +//! \f$(x-x_0)^2/B^2 + (\sqrt{(y-y_0)^2 + (z-z_0)^2} - A)^2/C^2 -1 \f$ +//============================================================================== + +class SurfaceXTorus : public CSGSurface { +public: + explicit SurfaceXTorus(pugi::xml_node surf_node); + double evaluate(Position r) const; + double distance(Position r, Direction u, bool coincident) const; + Direction normal(Position r) const; + void to_hdf5_inner(hid_t group_id) const; + + double x0_, y0_, z0_, A_, B_, C_; +}; + +//============================================================================== +//! A toroidal surface described by the quartic torus lies in the y direction +// +//! \f$(y-y_0)^2/B^2 + (\sqrt{(x-x_0)^2 + (z-z_0)^2} - A)^2/C^2 -1 \f$ +//============================================================================== + +class SurfaceYTorus : public CSGSurface { +public: + explicit SurfaceYTorus(pugi::xml_node surf_node); + double evaluate(Position r) const; + double distance(Position r, Direction u, bool coincident) const; + Direction normal(Position r) const; + void to_hdf5_inner(hid_t group_id) const; + + double x0_, y0_, z0_, A_, B_, C_; +}; + +//============================================================================== +//! A toroidal surface described by the quartic torus lies in the z direction +// +//! \f$(z-z_0)^2/B^2 + (\sqrt{(x-x_0)^2 + (y-y_0)^2} - A)^2/C^2 -1 \f$ +//============================================================================== + +class SurfaceZTorus : public CSGSurface { +public: + explicit SurfaceZTorus(pugi::xml_node surf_node); + double evaluate(Position r) const; + double distance(Position r, Direction u, bool coincident) const; + Direction normal(Position r) const; + void to_hdf5_inner(hid_t group_id) const; + + double x0_, y0_, z0_, A_, B_, C_; +}; + //============================================================================== // Non-member functions //============================================================================== diff --git a/openmc/surface.py b/openmc/surface.py index 41477fd1c9..4a174a23b0 100644 --- a/openmc/surface.py +++ b/openmc/surface.py @@ -648,7 +648,7 @@ class Plane(PlaneMixin, Surface): The 'C' parameter for the plane. Defaults to 0. d : float, optional The 'D' parameter for the plane. Defaults to 0. - boundary_type : {'transmission, 'vacuum', 'reflective', 'white'}, optional + boundary_type : {'transmission, 'vacuum', 'reflective', 'periodic', 'white'}, optional Boundary condition that defines the behavior for particles hitting the surface. Defaults to transmissive boundary condition where particles freely pass through the surface. @@ -668,7 +668,7 @@ class Plane(PlaneMixin, Surface): The 'C' parameter for the plane d : float The 'D' parameter for the plane - boundary_type : {'transmission, 'vacuum', 'reflective', 'white'} + boundary_type : {'transmission, 'vacuum', 'reflective', 'periodic', 'white'} Boundary condition that defines the behavior for particles hitting the surface. periodic_surface : openmc.Surface @@ -2064,7 +2064,7 @@ class Quadric(QuadricMixin, Surface): ---------- a, b, c, d, e, f, g, h, j, k : float, optional coefficients for the surface. All default to 0. - boundary_type : {'transmission, 'vacuum', 'reflective', 'periodic', 'white'}, optional + boundary_type : {'transmission, 'vacuum', 'reflective', 'white'}, optional Boundary condition that defines the behavior for particles hitting the surface. Defaults to transmissive boundary condition where particles freely pass through the surface. @@ -2078,7 +2078,7 @@ class Quadric(QuadricMixin, Surface): ---------- a, b, c, d, e, f, g, h, j, k : float coefficients for the surface - boundary_type : {'transmission, 'vacuum', 'reflective', 'periodic', 'white'} + boundary_type : {'transmission, 'vacuum', 'reflective', 'white'} Boundary condition that defines the behavior for particles hitting the surface. coefficients : dict @@ -2118,6 +2118,285 @@ class Quadric(QuadricMixin, Surface): return tuple(getattr(self, c) for c in self._coeff_keys) +class TorusMixin: + """A Mixin class implementing common functionality for torus surfaces""" + _coeff_keys = ('x0', 'y0', 'z0', 'a', 'b', 'c') + + def __init__(self, x0=0., y0=0., z0=0., a=0., b=0., c=0., **kwargs): + super().__init__(**kwargs) + for key, val in zip(self._coeff_keys, (x0, y0, z0, a, b, c)): + setattr(self, key, val) + + x0 = SurfaceCoefficient('x0') + y0 = SurfaceCoefficient('y0') + z0 = SurfaceCoefficient('z0') + a = SurfaceCoefficient('a') + b = SurfaceCoefficient('b') + c = SurfaceCoefficient('c') + + def translate(self, vector, inplace=False): + surf = self if inplace else self.clone() + surf.x0 += vector[0] + surf.y0 += vector[1] + surf.z0 += vector[2] + return surf + + def rotate(self, rotation, pivot=(0., 0., 0.), order='xyz', inplace=False): + pivot = np.asarray(pivot) + rotation = np.asarray(rotation, dtype=float) + + # Allow rotation matrix to be passed in directly, otherwise build it + if rotation.ndim == 2: + check_length('surface rotation', rotation.ravel(), 9) + Rmat = rotation + else: + Rmat = get_rotation_matrix(rotation, order=order) + + # Only can handle trivial rotation matrices + close = np.isclose + if not np.all(close(Rmat, -1.0) | close(Rmat, 0.0) | close(Rmat, 1.0)): + raise NotImplementedError('Torus surfaces cannot handle generic rotations') + + # Translate surface to pivot + surf = self.translate(-pivot, inplace=inplace) + + # Determine "center" of torus and a point above it (along main axis) + center = [surf.x0, surf.y0, surf.z0] + above_center = center.copy() + index = ['x-torus', 'y-torus', 'z-torus'].index(surf._type) + above_center[index] += 1 + + # Compute new rotated torus center + center = Rmat @ center + + # Figure out which axis should be used after rotation + above_center = Rmat @ above_center + new_index = np.where(np.isclose(np.abs(above_center - center), 1.0))[0][0] + cls = [XTorus, YTorus, ZTorus][new_index] + + # Create rotated torus + kwargs = { + 'boundary_type': surf.boundary_type, 'name': surf.name, + 'a': surf.a, 'b': surf.b, 'c': surf.c + } + if inplace: + kwargs['surface_id'] = surf.id + surf = cls(x0=center[0], y0=center[1], z0=center[2], **kwargs) + + return surf.translate(pivot, inplace=inplace) + + def _get_base_coeffs(self): + raise NotImplementedError + + +class XTorus(TorusMixin, Surface): + r"""A torus of the form :math:`(x - x_0)^2/B^2 + (\sqrt{(y - y_0)^2 + (z - + z_0)^2} - A)^2/C^2 - 1 = 0`. + + Parameters + ---------- + x0 : float + x-coordinate of the center of the axis of revolution + y0 : float + y-coordinate of the center of the axis of revolution + z0 : float + z-coordinate of the center of the axis of revolution + a : float + Major radius of the torus + b : float + Minor radius of the torus (parallel to axis of revolution) + c : float + Minor radius of the torus (perpendicular to axis of revolution) + kwargs : dict + Keyword arguments passed to the :class:`Surface` constructor + + Attributes + ---------- + x0 : float + x-coordinate of the center of the axis of revolution + y0 : float + y-coordinate of the center of the axis of revolution + z0 : float + z-coordinate of the center of the axis of revolution + a : float + Major radius of the torus + b : float + Minor radius of the torus (parallel to axis of revolution) + c : float + Minor radius of the torus (perpendicular to axis of revolution) + boundary_type : {'transmission, 'vacuum', 'reflective', 'white'} + Boundary condition that defines the behavior for particles hitting the + surface. + coefficients : dict + Dictionary of surface coefficients + id : int + Unique identifier for the surface + name : str + Name of the surface + type : str + Type of the surface + + """ + _type = 'x-torus' + + def evaluate(self, point): + x = point[0] - self.x0 + y = point[1] - self.y0 + z = point[2] - self.z0 + a = self.a + b = self.b + c = self.c + return (x*x)/(b*b) + (math.sqrt(y*y + z*z) - a)**2/(c*c) - 1 + + def bounding_box(self, side): + x0, y0, z0 = self.x0, self.y0, self.z0 + a, b, c = self.a, self.b, self.c + if side == '-': + return (np.array([x0 - b, y0 - a - c, z0 - a - c]), + np.array([x0 + b, y0 + a + c, z0 + a + c])) + elif side == '+': + return (np.array([-np.inf, -np.inf, -np.inf]), + np.array([np.inf, np.inf, np.inf])) + +class YTorus(TorusMixin, Surface): + r"""A torus of the form :math:`(y - y_0)^2/B^2 + (\sqrt{(x - x_0)^2 + (z - + z_0)^2} - A)^2/C^2 - 1 = 0`. + + Parameters + ---------- + x0 : float + x-coordinate of the center of the axis of revolution + y0 : float + y-coordinate of the center of the axis of revolution + z0 : float + z-coordinate of the center of the axis of revolution + a : float + Major radius of the torus + b : float + Minor radius of the torus (parallel to axis of revolution) + c : float + Minor radius of the torus (perpendicular to axis of revolution) + kwargs : dict + Keyword arguments passed to the :class:`Surface` constructor + + Attributes + ---------- + x0 : float + x-coordinate of the center of the axis of revolution + y0 : float + y-coordinate of the center of the axis of revolution + z0 : float + z-coordinate of the center of the axis of revolution + a : float + Major radius of the torus + b : float + Minor radius of the torus (parallel to axis of revolution) + c : float + Minor radius of the torus (perpendicular to axis of revolution) + boundary_type : {'transmission, 'vacuum', 'reflective', 'white'} + Boundary condition that defines the behavior for particles hitting the + surface. + coefficients : dict + Dictionary of surface coefficients + id : int + Unique identifier for the surface + name : str + Name of the surface + type : str + Type of the surface + + """ + _type = 'y-torus' + + def evaluate(self, point): + x = point[0] - self.x0 + y = point[1] - self.y0 + z = point[2] - self.z0 + a = self.a + b = self.b + c = self.c + return (y*y)/(b*b) + (math.sqrt(x*x + z*z) - a)**2/(c*c) - 1 + + def bounding_box(self, side): + x0, y0, z0 = self.x0, self.y0, self.z0 + a, b, c = self.a, self.b, self.c + if side == '-': + return (np.array([x0 - a - c, y0 - b, z0 - a - c]), + np.array([x0 + a + c, y0 + b, z0 + a + c])) + elif side == '+': + return (np.array([-np.inf, -np.inf, -np.inf]), + np.array([np.inf, np.inf, np.inf])) + +class ZTorus(TorusMixin, Surface): + r"""A torus of the form :math:`(z - z_0)^2/B^2 + (\sqrt{(x - x_0)^2 + (y - + y_0)^2} - A)^2/C^2 - 1 = 0`. + + Parameters + ---------- + x0 : float + x-coordinate of the center of the axis of revolution + y0 : float + y-coordinate of the center of the axis of revolution + z0 : float + z-coordinate of the center of the axis of revolution + a : float + Major radius of the torus + b : float + Minor radius of the torus (parallel to axis of revolution) + c : float + Minor radius of the torus (perpendicular to axis of revolution) + kwargs : dict + Keyword arguments passed to the :class:`Surface` constructor + + Attributes + ---------- + x0 : float + x-coordinate of the center of the axis of revolution + y0 : float + y-coordinate of the center of the axis of revolution + z0 : float + z-coordinate of the center of the axis of revolution + a : float + Major radius of the torus + b : float + Minor radius of the torus (parallel to axis of revolution) + c : float + Minor radius of the torus (perpendicular to axis of revolution) + boundary_type : {'transmission, 'vacuum', 'reflective', 'white'} + Boundary condition that defines the behavior for particles hitting the + surface. + coefficients : dict + Dictionary of surface coefficients + id : int + Unique identifier for the surface + name : str + Name of the surface + type : str + Type of the surface + """ + + _type = 'z-torus' + + def evaluate(self, point): + x = point[0] - self.x0 + y = point[1] - self.y0 + z = point[2] - self.z0 + a = self.a + b = self.b + c = self.c + return (z*z)/(b*b) + (math.sqrt(x*x + y*y) - a)**2/(c*c) - 1 + + def bounding_box(self, side): + x0, y0, z0 = self.x0, self.y0, self.z0 + a, b, c = self.a, self.b, self.c + if side == '-': + return (np.array([x0 - a - c, y0 - a - c, z0 - b]), + np.array([x0 + a + c, y0 + a + c, z0 + b])) + elif side == '+': + return (np.array([-np.inf, -np.inf, -np.inf]), + np.array([np.inf, np.inf, np.inf])) + + class Halfspace(Region): """A positive or negative half-space region. @@ -2360,7 +2639,7 @@ class Halfspace(Region): _SURFACE_CLASSES = {cls._type: cls for cls in Surface.__subclasses__()} -# Set virtual base classes for "casting" up the heirarchy +# Set virtual base classes for "casting" up the hierarchy Plane._virtual_base = Plane XPlane._virtual_base = Plane YPlane._virtual_base = Plane diff --git a/pyproject.toml b/pyproject.toml index d5970617a7..ac67355937 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -1,2 +1,2 @@ [build-system] -requires = ["setuptools", "wheel", "numpy", "cython"] +requires = ["setuptools", "wheel", "numpy<1.22", "cython"] diff --git a/src/external/LICENSE b/src/external/LICENSE new file mode 100644 index 0000000000..ac09c7fcb7 --- /dev/null +++ b/src/external/LICENSE @@ -0,0 +1,32 @@ +The quartic solver was obtained from the paper: Alberto Giacomo Orellana and +Cristiano De Michele, "Algorithm 1010: Boosting Efficiency in Solving Quartic +Equations with No Compromise in Accuracy," ACM Transactions on Mathematical +Software, 46 (2), pp. 1-28. https://doi.org/10.1145/3386241 + +OpenMC developers contacted the authors, who have agreed to license their +software under the simplified BSD license, reproduced below: + +------------------------------------------------------------------------------- +Copyright (c) 2020 Alberto Giacomo Orellana and Cristiano De Michele +All rights reserved. + +Redistribution and use in source and binary forms, with or without modification, +are permitted provided that the following conditions are met: + + 1. Redistributions of source code must retain the above copyright notice, this + list of conditions and the following disclaimer. + + 2. Redistributions in binary form must reproduce the above copyright notice, + this list of conditions and the following disclaimer in the documentation + and/or other materials provided with the distribution. + +THIS SOFTWARE IS PROVIDED BY THE PYNE DEVELOPMENT TEAM ``AS IS'' AND ANY EXPRESS +OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF +MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT +SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, +INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT +LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR +PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF +LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE +OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF +ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. diff --git a/src/external/quartic_solver.cpp b/src/external/quartic_solver.cpp new file mode 100644 index 0000000000..17410654fd --- /dev/null +++ b/src/external/quartic_solver.cpp @@ -0,0 +1,573 @@ +#include +#include +#include +#include + +namespace oqs { + +// pow(DBL_MAX,1.0/3.0)/1.618034; +constexpr double CUBIC_RESCAL_FACT = 3.488062113727083e+102; +// pow(DBL_MAX,1.0/4.0)/1.618034; +constexpr double QUART_RESCAL_FACT = 7.156344627944542e+76; +constexpr double MACHEPS = std::numeric_limits::epsilon(); + +double solve_cubic_analytic_depressed_handle_inf(double b, double c) +{ + /* find analytically the dominant root of a depressed cubic x^3+b*x+c + * where coefficients b and c are large (see sec. 2.2 in the manuscript) */ + double Q = -b / 3.0; + double R = 0.5 * c; + if (R == 0) { + return (b <= 0) ? std::sqrt(-b) : 0; + } + + double KK; + if (std::abs(Q) < std::abs(R)) { + double QR = Q / R; + double QRSQ = QR * QR; + KK = 1.0 - Q * QRSQ; + } else { + double RQ = R / Q; + KK = std::copysign(1.0, Q) * (RQ * RQ / Q - 1.0); + } + + if (KK < 0.0) { + double sqrtQ = std::sqrt(Q); + double theta = std::acos((R / std::abs(Q)) / sqrtQ); + if (2.0 * theta < M_PI) + return -2.0 * sqrtQ * std::cos(theta / 3.0); + else + return -2.0 * sqrtQ * std::cos((theta + 2.0 * M_PI) / 3.0); + } else { + double A; + if (std::abs(Q) < std::abs(R)) + A = -std::copysign(1.0, R) * cbrt(std::abs(R) * (1.0 + std::sqrt(KK))); + else { + A = -std::copysign(1.0, R) * + cbrt( + std::abs(R) + std::sqrt(std::abs(Q)) * std::abs(Q) * std::sqrt(KK)); + } + double B = (A == 0.0) ? 0.0 : Q / A; + return A + B; + } +} + +double solve_cubic_analytic_depressed(double b, double c) +{ + /* find analytically the dominant root of a depressed cubic x^3+b*x+c + * (see sec. 2.2 in the manuscript) */ + double Q = -b / 3.0; + double R = 0.5 * c; + if (std::abs(Q) > 1e102 || std::abs(R) > 1e154) { + return oqs::solve_cubic_analytic_depressed_handle_inf(b, c); + } + double Q3 = Q * Q * Q; + double R2 = R * R; + if (R2 < Q3) { + double theta = std::acos(R / std::sqrt(Q3)); + double sqrtQ = -2.0 * std::sqrt(Q); + if (2.0 * theta < M_PI) + return sqrtQ * std::cos(theta / 3.0); + else + return sqrtQ * std::cos((theta + 2.0 * M_PI) / 3.0); + } else { + double A = -std::copysign(1.0, R) * + std::pow(std::abs(R) + std::sqrt(R2 - Q3), 1.0 / 3.0); + double B = (A == 0.0) ? 0.0 : Q / A; + return A + B; /* this is always largest root even if A=B */ + } +} + +double calc_phi0(double a, double b, double c, double d, int scaled) +{ + /* find phi0 as the dominant root of the depressed and shifted cubic + * in eq. (79) (see also the discussion in sec. 2.2 of the manuscript) */ + double diskr = 9 * a * a - 24 * b; + /* eq. (87) */ + double s; + if (diskr > 0.0) { + diskr = std::sqrt(diskr); + s = -2 * b / (3 * a + std::copysign(diskr, a)); + } else { + s = -a / 4; + } + /* eqs. (83) */ + double aq = a + 4 * s; + double bq = b + 3 * s * (a + 2 * s); + double cq = c + s * (2 * b + s * (3 * a + 4 * s)); + double dq = d + s * (c + s * (b + s * (a + s))); + double gg = bq * bq / 9; + double hh = aq * cq; + + double g = hh - 4 * dq - 3 * gg; /* eq. (85) */ + double h = + (8 * dq + hh - 2 * gg) * bq / 3 - cq * cq - dq * aq * aq; /* eq. (86) */ + double rmax = oqs::solve_cubic_analytic_depressed(g, h); + if (std::isnan(rmax) || std::isinf(rmax)) { + rmax = oqs::solve_cubic_analytic_depressed_handle_inf(g, h); + if ((std::isnan(rmax) || std::isinf(rmax)) && scaled) { + // try harder: rescale also the depressed cubic if quartic has been + // already rescaled + double rfact = CUBIC_RESCAL_FACT; + double rfactsq = rfact * rfact; + double ggss = gg / rfactsq; + double hhss = hh / rfactsq; + double dqss = dq / rfactsq; + double aqs = aq / rfact; + double bqs = bq / rfact; + double cqs = cq / rfact; + ggss = bqs * bqs / 9.0; + hhss = aqs * cqs; + g = hhss - 4.0 * dqss - 3.0 * ggss; + h = (8.0 * dqss + hhss - 2.0 * ggss) * bqs / 3 - cqs * (cqs / rfact) - + (dq / rfact) * aqs * aqs; + rmax = oqs::solve_cubic_analytic_depressed(g, h); + if (std::isnan(rmax) || std::isinf(rmax)) { + rmax = oqs::solve_cubic_analytic_depressed_handle_inf(g, h); + } + rmax *= rfact; + } + } + /* Newton-Raphson used to refine phi0 (see end of sec. 2.2 in the manuscript) + */ + double x = rmax; + double xsq = x * x; + double xxx = x * xsq; + double gx = g * x; + double f = x * (xsq + g) + h; + double maxtt = std::max(std::abs(xxx), std::abs(gx)); + if (std::abs(h) > maxtt) + maxtt = std::abs(h); + + if (std::abs(f) > MACHEPS * maxtt) { + for (int iter = 0; iter < 8; iter++) { + double df = 3.0 * xsq + g; + if (df == 0) { + break; + } + double xold = x; + x += -f / df; + double fold = f; + xsq = x * x; + f = x * (xsq + g) + h; + if (f == 0) { + break; + } + + if (std::abs(f) >= std::abs(fold)) { + x = xold; + break; + } + } + } + return x; +} + +double calc_err_ldlt( + double b, double c, double d, double d2, double l1, double l2, double l3) +{ + /* Eqs. (29) and (30) in the manuscript */ + double sum = (b == 0) ? std::abs(d2 + l1 * l1 + 2.0 * l3) + : std::abs(((d2 + l1 * l1 + 2.0 * l3) - b) / b); + sum += (c == 0) ? std::abs(2.0 * d2 * l2 + 2.0 * l1 * l3) + : std::abs(((2.0 * d2 * l2 + 2.0 * l1 * l3) - c) / c); + sum += (d == 0) ? std::abs(d2 * l2 * l2 + l3 * l3) + : std::abs(((d2 * l2 * l2 + l3 * l3) - d) / d); + return sum; +} + +double calc_err_abcd_cmplx(double a, double b, double c, double d, + std::complex aq, std::complex bq, std::complex cq, + std::complex dq) +{ + /* Eqs. (68) and (69) in the manuscript for complex alpha1 (aq), beta1 (bq), + * alpha2 (cq) and beta2 (dq) */ + double sum = (d == 0) ? std::abs(bq * dq) : std::abs((bq * dq - d) / d); + sum += (c == 0) ? std::abs(bq * cq + aq * dq) + : std::abs(((bq * cq + aq * dq) - c) / c); + sum += (b == 0) ? std::abs(bq + aq * cq + dq) + : std::abs(((bq + aq * cq + dq) - b) / b); + sum += (a == 0) ? std::abs(aq + cq) : std::abs(((aq + cq) - a) / a); + return sum; +} +double calc_err_abcd(double a, double b, double c, double d, double aq, + double bq, double cq, double dq) +{ + /* Eqs. (68) and (69) in the manuscript for real alpha1 (aq), beta1 (bq), + * alpha2 (cq) and beta2 (dq)*/ + double sum = (d == 0) ? std::abs(bq * dq) : std::abs((bq * dq - d) / d); + sum += (c == 0) ? std::abs(bq * cq + aq * dq) + : std::abs(((bq * cq + aq * dq) - c) / c); + sum += (b == 0) ? std::abs(bq + aq * cq + dq) + : std::abs(((bq + aq * cq + dq) - b) / b); + sum += (a == 0) ? std::abs(aq + cq) : std::abs(((aq + cq) - a) / a); + return sum; +} + +double calc_err_abc( + double a, double b, double c, double aq, double bq, double cq, double dq) +{ + /* Eqs. (48)-(51) in the manuscript */ + double sum = (c == 0) ? std::abs(bq * cq + aq * dq) + : std::abs(((bq * cq + aq * dq) - c) / c); + sum += (b == 0) ? std::abs(bq + aq * cq + dq) + : std::abs(((bq + aq * cq + dq) - b) / b); + sum += (a == 0) ? std::abs(aq + cq) : std::abs(((aq + cq) - a) / a); + return sum; +} +void NRabcd(double a, double b, double c, double d, double* AQ, double* BQ, + double* CQ, double* DQ) +{ + /* Newton-Raphson described in sec. 2.3 of the manuscript for complex + * coefficients a,b,c,d */ + double xold[4], x[4], dx[4], det, Jinv[4][4], fvec[4], vr[4]; + x[0] = *AQ; + x[1] = *BQ; + x[2] = *CQ; + x[3] = *DQ; + vr[0] = d; + vr[1] = c; + vr[2] = b; + vr[3] = a; + fvec[0] = x[1] * x[3] - d; + fvec[1] = x[1] * x[2] + x[0] * x[3] - c; + fvec[2] = x[1] + x[0] * x[2] + x[3] - b; + fvec[3] = x[0] + x[2] - a; + double errf = 0; + for (int k1 = 0; k1 < 4; k1++) { + errf += (vr[k1] == 0) ? std::abs(fvec[k1]) : std::abs(fvec[k1] / vr[k1]); + } + for (int iter = 0; iter < 8; iter++) { + double x02 = x[0] - x[2]; + det = x[1] * x[1] + x[1] * (-x[2] * x02 - 2.0 * x[3]) + + x[3] * (x[0] * x02 + x[3]); + if (det == 0.0) + break; + Jinv[0][0] = x02; + Jinv[0][1] = x[3] - x[1]; + Jinv[0][2] = x[1] * x[2] - x[0] * x[3]; + Jinv[0][3] = -x[1] * Jinv[0][1] - x[0] * Jinv[0][2]; + Jinv[1][0] = x[0] * Jinv[0][0] + Jinv[0][1]; + Jinv[1][1] = -x[1] * Jinv[0][0]; + Jinv[1][2] = -x[1] * Jinv[0][1]; + Jinv[1][3] = -x[1] * Jinv[0][2]; + Jinv[2][0] = -Jinv[0][0]; + Jinv[2][1] = -Jinv[0][1]; + Jinv[2][2] = -Jinv[0][2]; + Jinv[2][3] = Jinv[0][2] * x[2] + Jinv[0][1] * x[3]; + Jinv[3][0] = -x[2] * Jinv[0][0] - Jinv[0][1]; + Jinv[3][1] = Jinv[0][0] * x[3]; + Jinv[3][2] = x[3] * Jinv[0][1]; + Jinv[3][3] = x[3] * Jinv[0][2]; + for (int k1 = 0; k1 < 4; k1++) { + dx[k1] = 0; + for (int k2 = 0; k2 < 4; k2++) + dx[k1] += Jinv[k1][k2] * fvec[k2]; + } + for (int k1 = 0; k1 < 4; k1++) + xold[k1] = x[k1]; + + for (int k1 = 0; k1 < 4; k1++) { + x[k1] += -dx[k1] / det; + } + fvec[0] = x[1] * x[3] - d; + fvec[1] = x[1] * x[2] + x[0] * x[3] - c; + fvec[2] = x[1] + x[0] * x[2] + x[3] - b; + fvec[3] = x[0] + x[2] - a; + double errfold = errf; + errf = 0; + for (int k1 = 0; k1 < 4; k1++) { + errf += (vr[k1] == 0) ? std::abs(fvec[k1]) : std::abs(fvec[k1] / vr[k1]); + } + if (errf == 0) + break; + if (errf >= errfold) { + for (int k1 = 0; k1 < 4; k1++) + x[k1] = xold[k1]; + break; + } + } + *AQ = x[0]; + *BQ = x[1]; + *CQ = x[2]; + *DQ = x[3]; +} + +void solve_quadratic(double a, double b, std::complex roots[2]) +{ + double diskr = a * a - 4 * b; + if (diskr >= 0.0) { + double div = -a - std::copysign(std::sqrt(diskr), a); + double zmax = div / 2; + double zmin = (zmax == 0.0) ? 0.0 : b / zmax; + + roots[0] = std::complex(zmax, 0.0); + roots[1] = std::complex(zmin, 0.0); + } else { + double sqrtd = std::sqrt(-diskr); + roots[0] = std::complex(-a / 2, sqrtd / 2); + roots[1] = std::complex(-a / 2, -sqrtd / 2); + } +} + +void quartic_solver(double coeff[5], std::complex roots[4]) +{ + /* USAGE: + * + * This routine calculates the roots of the quartic equation + * + * coeff[4]*x^4 + coeff[3]*x^3 + coeff[2]*x^2 + coeff[1]*x + coeff[0] = 0 + * + * if coeff[4] != 0 + * + * the four roots will be stored in the complex array roots[] + * + * */ + std::complex acx, bcx, ccx, dcx; + double l2m[12], d2m[12], res[12]; + double errv[3], aqv[3], cqv[3]; + int realcase[2]; + + double a = coeff[3] / coeff[4]; + double b = coeff[2] / coeff[4]; + double c = coeff[1] / coeff[4]; + double d = coeff[0] / coeff[4]; + double phi0 = oqs::calc_phi0(a, b, c, d, 0); + + // simple polynomial rescaling + double rfact = 1.0; + if (std::isnan(phi0) || std::isinf(phi0)) { + rfact = QUART_RESCAL_FACT; + a /= rfact; + double rfactsq = rfact * rfact; + b /= rfactsq; + c /= rfactsq * rfact; + d /= rfactsq * rfactsq; + phi0 = oqs::calc_phi0(a, b, c, d, 1); + } + double l1 = a / 2; /* eq. (16) */ + double l3 = b / 6 + phi0 / 2; /* eq. (18) */ + double del2 = c - a * l3; /* defined just after eq. (27) */ + int nsol = 0; + double bl311 = + 2. * b / 3. - phi0 - l1 * l1; /* This is d2 as defined in eq. (20)*/ + double dml3l3 = + d - l3 * l3; /* dml3l3 is d3 as defined in eq. (15) with d2=0 */ + + /* Three possible solutions for d2 and l2 (see eqs. (28) and discussion which + * follows) */ + if (bl311 != 0.0) { + d2m[nsol] = bl311; + l2m[nsol] = del2 / (2.0 * d2m[nsol]); + res[nsol] = oqs::calc_err_ldlt(b, c, d, d2m[nsol], l1, l2m[nsol], l3); + nsol++; + } + if (del2 != 0) { + l2m[nsol] = 2 * dml3l3 / del2; + if (l2m[nsol] != 0) { + d2m[nsol] = del2 / (2 * l2m[nsol]); + res[nsol] = oqs::calc_err_ldlt(b, c, d, d2m[nsol], l1, l2m[nsol], l3); + nsol++; + } + + d2m[nsol] = bl311; + l2m[nsol] = 2.0 * dml3l3 / del2; + res[nsol] = oqs::calc_err_ldlt(b, c, d, d2m[nsol], l1, l2m[nsol], l3); + nsol++; + } + + double l2, d2; + if (nsol == 0) { + l2 = d2 = 0.0; + } else { + /* we select the (d2,l2) pair which minimizes errors */ + double resmin; + int kmin; + for (int k1 = 0; k1 < nsol; k1++) { + if (k1 == 0 || res[k1] < resmin) { + resmin = res[k1]; + kmin = k1; + } + } + d2 = d2m[kmin]; + l2 = l2m[kmin]; + } + int whichcase = 0; + double aq, bq, cq, dq; + if (d2 < 0.0) { + /* Case I eqs. (37)-(40) */ + double gamma = std::sqrt(-d2); + aq = l1 + gamma; + bq = l3 + gamma * l2; + + cq = l1 - gamma; + dq = l3 - gamma * l2; + if (std::abs(dq) < std::abs(bq)) + dq = d / bq; + else if (std::abs(dq) > std::abs(bq)) + bq = d / dq; + if (std::abs(aq) < std::abs(cq)) { + nsol = 0; + if (dq != 0) { + aqv[nsol] = (c - bq * cq) / dq; /* see eqs. (47) */ + errv[nsol] = oqs::calc_err_abc(a, b, c, aqv[nsol], bq, cq, dq); + nsol++; + } + if (cq != 0) { + aqv[nsol] = (b - dq - bq) / cq; /* see eqs. (47) */ + errv[nsol] = oqs::calc_err_abc(a, b, c, aqv[nsol], bq, cq, dq); + nsol++; + } + aqv[nsol] = a - cq; /* see eqs. (47) */ + errv[nsol] = oqs::calc_err_abc(a, b, c, aqv[nsol], bq, cq, dq); + nsol++; + /* we select the value of aq (i.e. alpha1 in the manuscript) which + * minimizes errors */ + double errmin; + int kmin; + for (int k = 0; k < nsol; k++) { + if (k == 0 || errv[k] < errmin) { + kmin = k; + errmin = errv[k]; + } + } + aq = aqv[kmin]; + } else { + nsol = 0; + if (bq != 0) { + cqv[nsol] = (c - aq * dq) / bq; /* see eqs. (53) */ + errv[nsol] = oqs::calc_err_abc(a, b, c, aq, bq, cqv[nsol], dq); + nsol++; + } + if (aq != 0) { + cqv[nsol] = (b - bq - dq) / aq; /* see eqs. (53) */ + errv[nsol] = oqs::calc_err_abc(a, b, c, aq, bq, cqv[nsol], dq); + nsol++; + } + cqv[nsol] = a - aq; /* see eqs. (53) */ + errv[nsol] = oqs::calc_err_abc(a, b, c, aq, bq, cqv[nsol], dq); + nsol++; + /* we select the value of cq (i.e. alpha2 in the manuscript) which + * minimizes errors */ + double errmin; + int kmin; + for (int k = 0; k < nsol; k++) { + if (k == 0 || errv[k] < errmin) { + kmin = k; + errmin = errv[k]; + } + } + cq = cqv[kmin]; + } + realcase[0] = 1; + } else if (d2 > 0) { + /* Case II eqs. (53)-(56) */ + double gamma = std::sqrt(d2); + acx = std::complex(l1, gamma); + bcx = std::complex(l3, gamma * l2); + ccx = std::conj(acx); + dcx = std::conj(bcx); + realcase[0] = 0; + } else + realcase[0] = -1; // d2=0 + /* Case III: d2 is 0 or approximately 0 (in this case check which solution is + * better) */ + if (realcase[0] == -1 || + (std::abs(d2) <= + MACHEPS * (std::abs(2. * b / 3.) + std::abs(phi0) + l1 * l1))) { + double d3 = d - l3 * l3; + double err0 = 0.0; + if (realcase[0] == 1) + err0 = oqs::calc_err_abcd(a, b, c, d, aq, bq, cq, dq); + else if (realcase[0] == 0) + err0 = oqs::calc_err_abcd_cmplx(a, b, c, d, acx, bcx, ccx, dcx); + double aq1, bq1, cq1, dq1; + std::complex acx1, bcx1, ccx1, dcx1; + double err1 = 0.0; + if (d3 <= 0) { + realcase[1] = 1; + aq1 = l1; + bq1 = l3 + std::sqrt(-d3); + cq1 = l1; + dq1 = l3 - std::sqrt(-d3); + if (std::abs(dq1) < std::abs(bq1)) + dq1 = d / bq1; + else if (std::abs(dq1) > std::abs(bq1)) + bq1 = d / dq1; + err1 = oqs::calc_err_abcd(a, b, c, d, aq1, bq1, cq1, dq1); /* eq. (68) */ + } else { + /* complex */ + realcase[1] = 0; + acx1 = l1; + bcx1 = l3 + std::complex(0., std::sqrt(d3)); + ccx1 = l1; + dcx1 = std::conj(bcx1); + err1 = oqs::calc_err_abcd_cmplx(a, b, c, d, acx1, bcx1, ccx1, dcx1); + } + if (realcase[0] == -1 || err1 < err0) { + whichcase = 1; // d2 = 0 + if (realcase[1] == 1) { + aq = aq1; + bq = bq1; + cq = cq1; + dq = dq1; + } else { + acx = acx1; + bcx = bcx1; + ccx = ccx1; + dcx = dcx1; + } + } + } + if (realcase[whichcase] == 1) { + /* if alpha1, beta1, alpha2 and beta2 are real first refine + * the coefficient through a Newton-Raphson */ + oqs::NRabcd(a, b, c, d, &aq, &bq, &cq, &dq); + /* finally calculate the roots as roots of p1(x) and p2(x) (see end of + * sec. 2.1) */ + std::complex qroots[2]; + oqs::solve_quadratic(aq, bq, qroots); + roots[0] = qroots[0]; + roots[1] = qroots[1]; + oqs::solve_quadratic(cq, dq, qroots); + roots[2] = qroots[0]; + roots[3] = qroots[1]; + } else { + /* complex coefficients of p1 and p2 */ + if (whichcase == 0) { // d2!=0 + auto cdiskr = 0.25 * acx * acx - bcx; + /* calculate the roots as roots of p1(x) and p2(x) (see end of sec. 2.1) + */ + auto zx1 = -0.5 * acx + std::sqrt(cdiskr); + auto zx2 = -0.5 * acx - std::sqrt(cdiskr); + auto zxmax = (std::abs(zx1) > std::abs(zx2)) ? zx1 : zx2; + auto zxmin = bcx / zxmax; + roots[0] = zxmin; + roots[1] = std::conj(zxmin); + roots[2] = zxmax; + roots[3] = std::conj(zxmax); + } else { // d2 ~ 0 + /* never gets here! */ + auto cdiskr = std::sqrt(acx * acx - 4.0 * bcx); + auto zx1 = -0.5 * (acx + cdiskr); + auto zx2 = -0.5 * (acx - cdiskr); + auto zxmax = (std::abs(zx1) > std::abs(zx2)) ? zx1 : zx2; + auto zxmin = bcx / zxmax; + roots[0] = zxmax; + roots[1] = zxmin; + cdiskr = std::sqrt(ccx * ccx - 4.0 * dcx); + zx1 = -0.5 * (ccx + cdiskr); + zx2 = -0.5 * (ccx - cdiskr); + zxmax = (std::abs(zx1) > std::abs(zx2)) ? zx1 : zx2; + zxmin = dcx / zxmax; + roots[2] = zxmax; + roots[3] = zxmin; + } + } + if (rfact != 1.0) { + for (int k = 0; k < 4; k++) + roots[k] *= rfact; + } +} + +} // namespace oqs diff --git a/src/surface.cpp b/src/surface.cpp index 73a5aaaf08..dbfa0bf3dc 100644 --- a/src/surface.cpp +++ b/src/surface.cpp @@ -1,6 +1,7 @@ #include "openmc/surface.h" #include +#include #include #include @@ -10,6 +11,7 @@ #include "openmc/array.h" #include "openmc/container_util.h" #include "openmc/error.h" +#include "openmc/external/quartic_solver.h" #include "openmc/hdf5_interface.h" #include "openmc/math_functions.h" #include "openmc/random_lcg.h" @@ -45,7 +47,8 @@ void read_coeffs(pugi::xml_node surf_node, int surf_id, double& c1) // Parse the coefficients. int stat = sscanf(coeffs.c_str(), "%lf", &c1); if (stat != 1) { - fatal_error("Something went wrong reading surface coeffs"); + fatal_error(fmt::format( + "Something went wrong reading coeffs for surface {}", surf_id)); } } @@ -63,7 +66,8 @@ void read_coeffs( // Parse the coefficients. int stat = sscanf(coeffs.c_str(), "%lf %lf %lf", &c1, &c2, &c3); if (stat != 3) { - fatal_error("Something went wrong reading surface coeffs"); + fatal_error(fmt::format( + "Something went wrong reading coeffs for surface {}", surf_id)); } } @@ -81,7 +85,28 @@ void read_coeffs(pugi::xml_node surf_node, int surf_id, double& c1, double& c2, // Parse the coefficients. int stat = sscanf(coeffs.c_str(), "%lf %lf %lf %lf", &c1, &c2, &c3, &c4); if (stat != 4) { - fatal_error("Something went wrong reading surface coeffs"); + fatal_error(fmt::format( + "Something went wrong reading coeffs for surface {}", surf_id)); + } +} + +void read_coeffs(pugi::xml_node surf_node, int surf_id, double& c1, double& c2, + double& c3, double& c4, double& c5, double& c6) +{ + // Check the given number of coefficients. + std::string coeffs = get_node_value(surf_node, "coeffs"); + int n_words = word_count(coeffs); + if (n_words != 6) { + fatal_error(fmt::format( + "Surface {} expects 6 coeffs but was given {}", surf_id, n_words)); + } + + // Parse the coefficients. + int stat = sscanf( + coeffs.c_str(), "%lf %lf %lf %lf %lf %lf", &c1, &c2, &c3, &c4, &c5, &c6); + if (stat != 6) { + fatal_error(fmt::format( + "Something went wrong reading coeffs for surface {}", surf_id)); } } @@ -101,7 +126,8 @@ void read_coeffs(pugi::xml_node surf_node, int surf_id, double& c1, double& c2, int stat = sscanf(coeffs.c_str(), "%lf %lf %lf %lf %lf %lf %lf %lf %lf %lf", &c1, &c2, &c3, &c4, &c5, &c6, &c7, &c8, &c9, &c10); if (stat != 10) { - fatal_error("Something went wrong reading surface coeffs"); + fatal_error(fmt::format( + "Something went wrong reading coeffs for surface {}", surf_id)); } } @@ -983,6 +1009,210 @@ void SurfaceQuadric::to_hdf5_inner(hid_t group_id) const write_dataset(group_id, "coefficients", coeffs); } +//============================================================================== +// Torus helper functions +//============================================================================== + +double torus_distance(double x1, double x2, double x3, double u1, double u2, + double u3, double A, double B, double C, bool coincident) +{ + // Coefficients for equation: (c2 t^2 + c1 t + c0)^2 = c2' t^2 + c1' t + c0' + double D = (C * C) / (B * B); + double c2 = u1 * u1 + u2 * u2 + D * u3 * u3; + double c1 = 2 * (u1 * x1 + u2 * x2 + D * u3 * x3); + double c0 = x1 * x1 + x2 * x2 + D * x3 * x3 + A * A - C * C; + double four_A2 = 4 * A * A; + double c2p = four_A2 * (u1 * u1 + u2 * u2); + double c1p = 2 * four_A2 * (u1 * x1 + u2 * x2); + double c0p = four_A2 * (x1 * x1 + x2 * x2); + + // Coefficient for equation: a t^4 + b t^3 + c t^2 + d t + e = 0 + double coeff[5]; + coeff[0] = c0 * c0 - c0p; + coeff[1] = 2 * c0 * c1 - c1p; + coeff[2] = c1 * c1 + 2 * c0 * c2 - c2p; + coeff[3] = 2 * c1 * c2; + coeff[4] = c2 * c2; + + std::complex roots[4]; + oqs::quartic_solver(coeff, roots); + + // Find smallest positive, real root. In the case where the particle is + // coincident with the surface, we are sure to have one root very close to + // zero but possibly small and positive. A tolerance is set to discard that + // zero. + double distance = INFTY; + double cutoff = coincident ? 1e-10 : 0.0; + for (int i = 0; i < 4; ++i) { + if (roots[i].imag() == 0) { + double root = roots[i].real(); + if (root > cutoff && root < distance) { + distance = root; + } + } + } + return distance; +} + +//============================================================================== +// SurfaceXTorus implementation +//============================================================================== + +SurfaceXTorus::SurfaceXTorus(pugi::xml_node surf_node) : CSGSurface(surf_node) +{ + read_coeffs(surf_node, id_, x0_, y0_, z0_, A_, B_, C_); +} + +void SurfaceXTorus::to_hdf5_inner(hid_t group_id) const +{ + write_string(group_id, "type", "x-torus", false); + std::array coeffs {{x0_, y0_, z0_, A_, B_, C_}}; + write_dataset(group_id, "coefficients", coeffs); +} + +double SurfaceXTorus::evaluate(Position r) const +{ + double x = r.x - x0_; + double y = r.y - y0_; + double z = r.z - z0_; + return (x * x) / (B_ * B_) + + std::pow(std::sqrt(y * y + z * z) - A_, 2) / (C_ * C_) - 1.; +} + +double SurfaceXTorus::distance(Position r, Direction u, bool coincident) const +{ + double x = r.x - x0_; + double y = r.y - y0_; + double z = r.z - z0_; + return torus_distance(y, z, x, u.y, u.z, u.x, A_, B_, C_, coincident); +} + +Direction SurfaceXTorus::normal(Position r) const +{ + // reduce the expansion of the full form for torus + double x = r.x - x0_; + double y = r.y - y0_; + double z = r.z - z0_; + + // f(x,y,z) = x^2/B^2 + (sqrt(y^2 + z^2) - A)^2/C^2 - 1 + // ∂f/∂x = 2x/B^2 + // ∂f/∂y = 2y(g - A)/(g*C^2) where g = sqrt(y^2 + z^2) + // ∂f/∂z = 2z(g - A)/(g*C^2) + // Multiplying by g*C^2*B^2 / 2 gives: + double g = std::sqrt(y * y + z * z); + double nx = C_ * C_ * g * x; + double ny = y * (g - A_) * B_ * B_; + double nz = z * (g - A_) * B_ * B_; + Direction n(nx, ny, nz); + return n / n.norm(); +} + +//============================================================================== +// SurfaceYTorus implementation +//============================================================================== + +SurfaceYTorus::SurfaceYTorus(pugi::xml_node surf_node) : CSGSurface(surf_node) +{ + read_coeffs(surf_node, id_, x0_, y0_, z0_, A_, B_, C_); +} + +void SurfaceYTorus::to_hdf5_inner(hid_t group_id) const +{ + write_string(group_id, "type", "y-torus", false); + std::array coeffs {{x0_, y0_, z0_, A_, B_, C_}}; + write_dataset(group_id, "coefficients", coeffs); +} + +double SurfaceYTorus::evaluate(Position r) const +{ + double x = r.x - x0_; + double y = r.y - y0_; + double z = r.z - z0_; + return (y * y) / (B_ * B_) + + std::pow(std::sqrt(x * x + z * z) - A_, 2) / (C_ * C_) - 1.; +} + +double SurfaceYTorus::distance(Position r, Direction u, bool coincident) const +{ + double x = r.x - x0_; + double y = r.y - y0_; + double z = r.z - z0_; + return torus_distance(x, z, y, u.x, u.z, u.y, A_, B_, C_, coincident); +} + +Direction SurfaceYTorus::normal(Position r) const +{ + // reduce the expansion of the full form for torus + double x = r.x - x0_; + double y = r.y - y0_; + double z = r.z - z0_; + + // f(x,y,z) = y^2/B^2 + (sqrt(x^2 + z^2) - A)^2/C^2 - 1 + // ∂f/∂x = 2x(g - A)/(g*C^2) where g = sqrt(x^2 + z^2) + // ∂f/∂y = 2y/B^2 + // ∂f/∂z = 2z(g - A)/(g*C^2) + // Multiplying by g*C^2*B^2 / 2 gives: + double g = std::sqrt(x * x + z * z); + double nx = x * (g - A_) * B_ * B_; + double ny = C_ * C_ * g * y; + double nz = z * (g - A_) * B_ * B_; + Direction n(nx, ny, nz); + return n / n.norm(); +} + +//============================================================================== +// SurfaceZTorus implementation +//============================================================================== + +SurfaceZTorus::SurfaceZTorus(pugi::xml_node surf_node) : CSGSurface(surf_node) +{ + read_coeffs(surf_node, id_, x0_, y0_, z0_, A_, B_, C_); +} + +void SurfaceZTorus::to_hdf5_inner(hid_t group_id) const +{ + write_string(group_id, "type", "z-torus", false); + std::array coeffs {{x0_, y0_, z0_, A_, B_, C_}}; + write_dataset(group_id, "coefficients", coeffs); +} + +double SurfaceZTorus::evaluate(Position r) const +{ + double x = r.x - x0_; + double y = r.y - y0_; + double z = r.z - z0_; + return (z * z) / (B_ * B_) + + std::pow(std::sqrt(x * x + y * y) - A_, 2) / (C_ * C_) - 1.; +} + +double SurfaceZTorus::distance(Position r, Direction u, bool coincident) const +{ + double x = r.x - x0_; + double y = r.y - y0_; + double z = r.z - z0_; + return torus_distance(x, y, z, u.x, u.y, u.z, A_, B_, C_, coincident); +} + +Direction SurfaceZTorus::normal(Position r) const +{ + // reduce the expansion of the full form for torus + double x = r.x - x0_; + double y = r.y - y0_; + double z = r.z - z0_; + + // f(x,y,z) = z^2/B^2 + (sqrt(x^2 + y^2) - A)^2/C^2 - 1 + // ∂f/∂x = 2x(g - A)/(g*C^2) where g = sqrt(x^2 + y^2) + // ∂f/∂y = 2y(g - A)/(g*C^2) + // ∂f/∂z = 2z/B^2 + // Multiplying by g*C^2*B^2 / 2 gives: + double g = std::sqrt(x * x + y * y); + double nx = x * (g - A_) * B_ * B_; + double ny = y * (g - A_) * B_ * B_; + double nz = C_ * C_ * g * z; + Position n(nx, ny, nz); + return n / n.norm(); +} + //============================================================================== void read_surfaces(pugi::xml_node node) @@ -1042,6 +1272,15 @@ void read_surfaces(pugi::xml_node node) } else if (surf_type == "quadric") { model::surfaces.push_back(make_unique(surf_node)); + } else if (surf_type == "x-torus") { + model::surfaces.push_back(std::make_unique(surf_node)); + + } else if (surf_type == "y-torus") { + model::surfaces.push_back(std::make_unique(surf_node)); + + } else if (surf_type == "z-torus") { + model::surfaces.push_back(std::make_unique(surf_node)); + } else { fatal_error(fmt::format("Invalid surface type, \"{}\"", surf_type)); } diff --git a/tests/regression_tests/torus/__init__.py b/tests/regression_tests/torus/__init__.py new file mode 100644 index 0000000000..e69de29bb2 diff --git a/tests/regression_tests/torus/inputs_true.dat b/tests/regression_tests/torus/inputs_true.dat new file mode 100644 index 0000000000..0045d7e466 --- /dev/null +++ b/tests/regression_tests/torus/inputs_true.dat @@ -0,0 +1,34 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + eigenvalue + 1000 + 10 + 5 + diff --git a/tests/regression_tests/torus/results_true.dat b/tests/regression_tests/torus/results_true.dat new file mode 100644 index 0000000000..42444ad09a --- /dev/null +++ b/tests/regression_tests/torus/results_true.dat @@ -0,0 +1,2 @@ +k-combined: +7.628424E-01 2.557300E-02 diff --git a/tests/regression_tests/torus/test.py b/tests/regression_tests/torus/test.py new file mode 100644 index 0000000000..8970a5cc77 --- /dev/null +++ b/tests/regression_tests/torus/test.py @@ -0,0 +1,39 @@ +import openmc +import pytest + +from tests.testing_harness import PyAPITestHarness + + +@pytest.fixture +def model(): + model = openmc.Model() + fuel = openmc.Material() + fuel.set_density('g/cm3', 12.0) + fuel.add_nuclide('U235', 1.0) + al = openmc.Material() + al.set_density('g/cm3', 1.0) + al.add_nuclide('H1', 1.0) + model.materials.extend([fuel, al]) + + # 🍩🍩🍩 + zt = openmc.ZTorus(a=3, b=1.5, c=1) + xt = openmc.XTorus(x0=6, a=3, b=1.5, c=1) + yt = openmc.YTorus(x0=6, a=6, b=1, c=0.75) + box = openmc.model.RectangularParallelepiped(-5, 14, -5, 5, -8, 8, + boundary_type='vacuum') + + xt_cell = openmc.Cell(fill=fuel, region=-xt) + yt_cell = openmc.Cell(fill=fuel, region=-yt) + zt_cell = openmc.Cell(fill=fuel, region=-zt) + outer_cell = openmc.Cell(fill=al, region=-box & +xt & +yt & +zt) + model.geometry = openmc.Geometry([xt_cell, yt_cell, zt_cell, outer_cell]) + + model.settings.particles = 1000 + model.settings.batches = 10 + model.settings.inactive = 5 + return model + + +def test_torus(model): + harness = PyAPITestHarness('statepoint.10.h5', model) + harness.main() diff --git a/tests/unit_tests/test_surface.py b/tests/unit_tests/test_surface.py index e7d37c4e91..12cd8c9d20 100644 --- a/tests/unit_tests/test_surface.py +++ b/tests/unit_tests/test_surface.py @@ -605,3 +605,151 @@ def test_cylinder_from_points_axis(): assert (d, e, f) == pytest.approx((0., 0., 0.)) assert (g, h, j) == pytest.approx((-4., 0., -10.)) assert k == pytest.approx(13.) + + +def torus_common(center, R, r1, r2, cls): + x, y, z = center + s = cls(x0=x, y0=y, z0=z, a=R, b=r1, c=r2) + assert s.x0 == x + assert s.y0 == y + assert s.z0 == z + assert s.a == R + assert s.b == r1 + assert s.c == r2 + + # evaluate method + assert s.evaluate((x, y, z)) > 0.0 + + # translate method + trans = np.array([1.0, 1.5, -2.0]) + st = s.translate(trans) + assert st.x0 == s.x0 + trans[0] + assert st.y0 == s.y0 + trans[1] + assert st.z0 == s.z0 + trans[2] + assert st.a == s.a + assert st.b == s.b + assert st.c == s.c + + # trivial rotations + for rotation in [(0., 0., 0.), (180., 0., 0.), (0., 180., 0.), (0., 0., 180.)]: + sr = s.rotate(rotation) + assert type(sr) == type(s) + assert (sr.a, sr.b, sr.c) == (s.a, s.b, s.c) + + # can't do generic rotate at present + with pytest.raises(NotImplementedError): + s.rotate((0., 45., 0.)) + + # Check bounding box + ll, ur = (+s).bounding_box + assert np.all(np.isinf(ll)) + assert np.all(np.isinf(ur)) + + ll, ur = (-s).bounding_box + llt, urt = (-st).bounding_box + np.testing.assert_allclose(ll + trans, llt) + np.testing.assert_allclose(ur + trans, urt) + + # Make sure repr works + repr(s) + + return s + + +def test_xtorus(): + x, y, z = 2, -4, 5 + R, r1, r2 = 3, 1.5, 1 + s = torus_common((x, y, z), R, r1, r2, openmc.XTorus) + + # evaluate method (points inside torus) + assert s.evaluate((x, y + R, z)) < 0.0 + assert s.evaluate((x, y - R, z)) < 0.0 + assert s.evaluate((x, y, z + R)) < 0.0 + assert s.evaluate((x, y, z - R)) < 0.0 + + # evaluate method (points on torus) + assert s.evaluate((x, y + R + r2, z)) == pytest.approx(0.0) + assert s.evaluate((x, y - R - r2, z)) == pytest.approx(0.0) + assert s.evaluate((x, y, z + R - r2)) == pytest.approx(0.0) + assert s.evaluate((x, y, z - R + r2)) == pytest.approx(0.0) + assert s.evaluate((x + r1, y + R, z)) == pytest.approx(0.0) + + # evaluate method (points outside torus) + assert s.evaluate((x, y + R, z + R)) > 0.0 + assert s.evaluate((x, y - R, z + R)) > 0.0 + assert s.evaluate((x, y + R + r2 + 0.01, z)) > 0.0 + assert s.evaluate((x, y, z + R + r2 + 0.01)) > 0.0 + assert s.evaluate((x + r1 + 0.01, y, z + R)) > 0.0 + assert s.evaluate((x + r1 + 0.01, y + R, z)) > 0.0 + + # rotation + sr = s.rotate((0., 0., 90.)) + assert isinstance(sr, openmc.YTorus) + sr = s.rotate((0., 90., 0.)) + assert isinstance(sr, openmc.ZTorus) + + +def test_ytorus(): + x, y, z = 2, -4, 5 + R, r1, r2 = 3, 1.5, 1 + s = torus_common((x, y, z), R, r1, r2, openmc.YTorus) + + # evaluate method (points inside torus) + assert s.evaluate((x + R, y, z)) < 0.0 + assert s.evaluate((x - R, y, z)) < 0.0 + assert s.evaluate((x, y, z + R)) < 0.0 + assert s.evaluate((x, y, z - R)) < 0.0 + + # evaluate method (points on torus) + assert s.evaluate((x + R + r2, y, z)) == pytest.approx(0.0) + assert s.evaluate((x - R - r2, y, z)) == pytest.approx(0.0) + assert s.evaluate((x, y, z + R - r2)) == pytest.approx(0.0) + assert s.evaluate((x, y, z - R + r2)) == pytest.approx(0.0) + assert s.evaluate((x + R, y + r1, z)) == pytest.approx(0.0) + + # evaluate method (points outside torus) + assert s.evaluate((x + R, y, z + R)) > 0.0 + assert s.evaluate((x - R, y, z + R)) > 0.0 + assert s.evaluate((x + R + r2 + 0.01, y, z)) > 0.0 + assert s.evaluate((x, y, z + R + r2 + 0.01)) > 0.0 + assert s.evaluate((x, y + r1 + 0.01, z + R)) > 0.0 + assert s.evaluate((x + R, y + r1 + 0.01, z)) > 0.0 + + # rotation + sr = s.rotate((90., 0., 0.)) + assert isinstance(sr, openmc.ZTorus) + sr = s.rotate((0., 0., 90.)) + assert isinstance(sr, openmc.XTorus) + + +def test_ztorus(): + x, y, z = 2, -4, 5 + R, r1, r2 = 3, 1.5, 1 + s = torus_common((x, y, z), R, r1, r2, openmc.ZTorus) + + # evaluate method (points inside torus) + assert s.evaluate((x, y + R, z)) < 0.0 + assert s.evaluate((x, y - R, z)) < 0.0 + assert s.evaluate((x + R, y, z)) < 0.0 + assert s.evaluate((x - R, y, z)) < 0.0 + + # evaluate method (points on torus) + assert s.evaluate((x, y + R + r2, z)) == pytest.approx(0.0) + assert s.evaluate((x, y - R - r2, z)) == pytest.approx(0.0) + assert s.evaluate((x + R - r2, y, z)) == pytest.approx(0.0) + assert s.evaluate((x - R + r2, y, z)) == pytest.approx(0.0) + assert s.evaluate((x, y + R, z + r1)) == pytest.approx(0.0) + + # evaluate method (points outside torus) + assert s.evaluate((x + R, y + R, z)) > 0.0 + assert s.evaluate((x + R, y - R, z)) > 0.0 + assert s.evaluate((x, y + R + r2 + 0.01, z)) > 0.0 + assert s.evaluate((x + R + r2 + 0.01, y, z)) > 0.0 + assert s.evaluate((x + R, y, z + r1 + 0.01)) > 0.0 + assert s.evaluate((x, y + R, z + r1 + 0.01)) > 0.0 + + # rotation + sr = s.rotate((90., 0., 0.)) + assert isinstance(sr, openmc.YTorus) + sr = s.rotate((0., 90., 0.)) + assert isinstance(sr, openmc.XTorus) diff --git a/tests/unit_tests/test_torus.py b/tests/unit_tests/test_torus.py new file mode 100644 index 0000000000..a5f49faae2 --- /dev/null +++ b/tests/unit_tests/test_torus.py @@ -0,0 +1,47 @@ +from itertools import combinations +from random import uniform +import openmc +import pytest + + +def get_torus_keff(cls, R, r, center=(0, 0, 0)): + model = openmc.Model() + mat = openmc.Material() + mat.add_nuclide('U235', 1.0) + mat.set_density('g/cm3', 10.0) + + x, y, z = center + torus = cls(x0=x, y0=y, z0=z, a=R, b=r, c=r) + sphere = openmc.Sphere(x0=x, y0=y, z0=z, r=R + r + 1, boundary_type="vacuum") + torus_cell = openmc.Cell(fill=mat, region=-torus) + outer_cell = openmc.Cell(region=+torus & -sphere) + model.geometry = openmc.Geometry([torus_cell, outer_cell]) + + model.settings.source = openmc.Source(space=openmc.stats.Point(center)) + model.settings.batches = 10 + model.settings.inactive = 5 + model.settings.particles = 1000 + + sp_path = model.run() + with openmc.StatePoint(sp_path) as sp: + return sp.k_combined + + +@pytest.mark.parametrize("R,r", [(2.1, 2.0), (3.0, 1.0)]) +def test_torus_keff(R, r, run_in_tmpdir): + random_point = lambda: (uniform(-5, 5), uniform(-5, 5), uniform(-5, 5)) + keffs = [ + get_torus_keff(openmc.XTorus, R, r), + get_torus_keff(openmc.XTorus, R, r, random_point()), + get_torus_keff(openmc.YTorus, R, r), + get_torus_keff(openmc.YTorus, R, r, random_point()), + get_torus_keff(openmc.ZTorus, R, r), + get_torus_keff(openmc.ZTorus, R, r, random_point()) + ] + + # For each combination of keff values, their difference should be within + # uncertainty (3 std dev) + for k1, k2 in combinations(keffs, 2): + print(k1, k2) + diff = k1 - k2 + assert abs(diff.n) < 3*diff.s diff --git a/tools/ci/gha-install.sh b/tools/ci/gha-install.sh index 8c4762e3e4..aa40eb90b1 100755 --- a/tools/ci/gha-install.sh +++ b/tools/ci/gha-install.sh @@ -1,10 +1,13 @@ #!/bin/bash set -ex -# Upgrade pip, pytest, numpy before doing anything else +# Upgrade pip, pytest, numpy before doing anything else. +# TODO: numpy 1.22 results in several failing tests, so we force a lower version +# for now (similar change made in pyproject.toml). When this is removed, those +# tests will need to be updated. pip install --upgrade pip pip install --upgrade pytest -pip install --upgrade numpy +pip install --upgrade "numpy<1.22" # Install NJOY 2016 ./tools/ci/gha-install-njoy.sh