Add documentation on torus distance/normal

docs/source/io_formats/geometry.rst

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

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.

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

src/surface.cpp

@@ -1013,21 +1013,21 @@ void SurfaceQuadric::to_hdf5_inner(hid_t group_id) const
 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 = d2 t^2 + d1 t + d0
+  // 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 d2 = four_A2 * (u1 * u1 + u2 * u2);
-  double d1 = 2 * four_A2 * (u1 * x1 + u2 * x2);
-  double d0 = four_A2 * (x1 * x1 + x2 * x2);
+  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 - d0;
-  coeff[1] = 2 * c0 * c1 - d1;
-  coeff[2] = c1 * c1 + 2 * c0 * c2 - d2;
+  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;