mirror of
https://github.com/openmc-dev/openmc.git
synced 2026-07-28 06:05:58 -04:00
Merge pull request #1933 from paulromano/torus-oqs
Implementation of torus surfaces
This commit is contained in:
Patrick Shriwise
GitHub
commit
0694afc554
22 changed files with 1603 additions and 27 deletions
|
|
@ -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)
|
||||
|
|
|
|||
|
|
@ -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.
|
||||
|
|
|
|||
12
README.md
12
README.md
|
|
@ -1,6 +1,6 @@
|
|||
# OpenMC Monte Carlo Particle Transport Code
|
||||
|
||||
[](http://openmc.readthedocs.io/en/latest/license.html)
|
||||
[](https://docs.openmc.org/en/latest/license.html)
|
||||
[](https://github.com/openmc-dev/openmc/actions?query=workflow%3ACI)
|
||||
[](https://coveralls.io/github/openmc-dev/openmc?branch=develop)
|
||||
[](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).
|
||||
|
|
|
|||
|
|
@ -26,7 +26,7 @@ Each ``<surface>`` 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
|
||||
|
||||
|
|
|
|||
|
|
@ -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.
|
||||
|
|
|
|||
|
|
@ -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 <https://doi.org/10.1145/3386241>`_ 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:
|
||||
|
||||
|
|
|
|||
|
|
@ -80,6 +80,9 @@ Building geometry
|
|||
openmc.YCone
|
||||
openmc.ZCone
|
||||
openmc.Quadric
|
||||
openmc.XTorus
|
||||
openmc.YTorus
|
||||
openmc.ZTorus
|
||||
openmc.Halfspace
|
||||
openmc.Intersection
|
||||
openmc.Union
|
||||
|
|
|
|||
|
|
@ -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
|
||||
|
|
|
|||
10
include/openmc/external/quartic_solver.h
vendored
Normal file
10
include/openmc/external/quartic_solver.h
vendored
Normal file
|
|
@ -0,0 +1,10 @@
|
|||
#ifndef OPENMC_EXTERNAL_QUARTIC_SOLVER_H
|
||||
#define OPENMC_EXTERNAL_QUARTIC_SOLVER_H
|
||||
|
||||
#include <complex>
|
||||
|
||||
namespace oqs {
|
||||
void quartic_solver(double coeff[5], std::complex<double> roots[4]);
|
||||
} // end namespace oqs
|
||||
|
||||
#endif // OPENMC_EXTERNAL_QUARTIC_SOLVER_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
|
||||
//==============================================================================
|
||||
|
|
|
|||
|
|
@ -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
|
||||
|
|
|
|||
|
|
@ -1,2 +1,2 @@
|
|||
[build-system]
|
||||
requires = ["setuptools", "wheel", "numpy", "cython"]
|
||||
requires = ["setuptools", "wheel", "numpy<1.22", "cython"]
|
||||
|
|
|
|||
32
src/external/LICENSE
vendored
Normal file
32
src/external/LICENSE
vendored
Normal file
|
|
@ -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.
|
||||
573
src/external/quartic_solver.cpp
vendored
Normal file
573
src/external/quartic_solver.cpp
vendored
Normal file
|
|
@ -0,0 +1,573 @@
|
|||
#include <algorithm>
|
||||
#include <cmath>
|
||||
#include <complex>
|
||||
#include <cstdlib>
|
||||
|
||||
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<double>::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<double> aq, std::complex<double> bq, std::complex<double> cq,
|
||||
std::complex<double> 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<double> 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<double>(zmax, 0.0);
|
||||
roots[1] = std::complex<double>(zmin, 0.0);
|
||||
} else {
|
||||
double sqrtd = std::sqrt(-diskr);
|
||||
roots[0] = std::complex<double>(-a / 2, sqrtd / 2);
|
||||
roots[1] = std::complex<double>(-a / 2, -sqrtd / 2);
|
||||
}
|
||||
}
|
||||
|
||||
void quartic_solver(double coeff[5], std::complex<double> 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<double> 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<double>(l1, gamma);
|
||||
bcx = std::complex<double>(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<double> 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<double>(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<double> 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
|
||||
247
src/surface.cpp
247
src/surface.cpp
|
|
@ -1,6 +1,7 @@
|
|||
#include "openmc/surface.h"
|
||||
|
||||
#include <cmath>
|
||||
#include <complex>
|
||||
#include <set>
|
||||
#include <utility>
|
||||
|
||||
|
|
@ -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<double> 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<double, 6> 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<double, 6> 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<double, 6> 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<SurfaceQuadric>(surf_node));
|
||||
|
||||
} else if (surf_type == "x-torus") {
|
||||
model::surfaces.push_back(std::make_unique<SurfaceXTorus>(surf_node));
|
||||
|
||||
} else if (surf_type == "y-torus") {
|
||||
model::surfaces.push_back(std::make_unique<SurfaceYTorus>(surf_node));
|
||||
|
||||
} else if (surf_type == "z-torus") {
|
||||
model::surfaces.push_back(std::make_unique<SurfaceZTorus>(surf_node));
|
||||
|
||||
} else {
|
||||
fatal_error(fmt::format("Invalid surface type, \"{}\"", surf_type));
|
||||
}
|
||||
|
|
|
|||
0
tests/regression_tests/torus/__init__.py
Normal file
0
tests/regression_tests/torus/__init__.py
Normal file
34
tests/regression_tests/torus/inputs_true.dat
Normal file
34
tests/regression_tests/torus/inputs_true.dat
Normal file
|
|
@ -0,0 +1,34 @@
|
|||
<?xml version='1.0' encoding='utf-8'?>
|
||||
<geometry>
|
||||
<cell id="1" material="1" region="-2" universe="1" />
|
||||
<cell id="2" material="1" region="-3" universe="1" />
|
||||
<cell id="3" material="1" region="-1" universe="1" />
|
||||
<cell id="4" material="2" region="4 -5 6 -7 8 -9 2 3 1" universe="1" />
|
||||
<surface coeffs="0.0 0.0 0.0 3 1.5 1" id="1" type="z-torus" />
|
||||
<surface coeffs="6 0.0 0.0 3 1.5 1" id="2" type="x-torus" />
|
||||
<surface coeffs="6 0.0 0.0 6 1 0.75" id="3" type="y-torus" />
|
||||
<surface boundary="vacuum" coeffs="-5" id="4" type="x-plane" />
|
||||
<surface boundary="vacuum" coeffs="14" id="5" type="x-plane" />
|
||||
<surface boundary="vacuum" coeffs="-5" id="6" type="y-plane" />
|
||||
<surface boundary="vacuum" coeffs="5" id="7" type="y-plane" />
|
||||
<surface boundary="vacuum" coeffs="-8" id="8" type="z-plane" />
|
||||
<surface boundary="vacuum" coeffs="8" id="9" type="z-plane" />
|
||||
</geometry>
|
||||
<?xml version='1.0' encoding='utf-8'?>
|
||||
<materials>
|
||||
<material depletable="true" id="1">
|
||||
<density units="g/cm3" value="12.0" />
|
||||
<nuclide ao="1.0" name="U235" />
|
||||
</material>
|
||||
<material id="2">
|
||||
<density units="g/cm3" value="1.0" />
|
||||
<nuclide ao="1.0" name="H1" />
|
||||
</material>
|
||||
</materials>
|
||||
<?xml version='1.0' encoding='utf-8'?>
|
||||
<settings>
|
||||
<run_mode>eigenvalue</run_mode>
|
||||
<particles>1000</particles>
|
||||
<batches>10</batches>
|
||||
<inactive>5</inactive>
|
||||
</settings>
|
||||
2
tests/regression_tests/torus/results_true.dat
Normal file
2
tests/regression_tests/torus/results_true.dat
Normal file
|
|
@ -0,0 +1,2 @@
|
|||
k-combined:
|
||||
7.628424E-01 2.557300E-02
|
||||
39
tests/regression_tests/torus/test.py
Normal file
39
tests/regression_tests/torus/test.py
Normal file
|
|
@ -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()
|
||||
|
|
@ -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)
|
||||
|
|
|
|||
47
tests/unit_tests/test_torus.py
Normal file
47
tests/unit_tests/test_torus.py
Normal file
|
|
@ -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
|
||||
|
|
@ -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
|
||||
|
|
|
|||
Loading…
Add table
Add a link
Reference in a new issue