diff --git a/docs/source/io_formats/geometry.rst b/docs/source/io_formats/geometry.rst index b04edca2b..ff2f6eb74 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 9b92b3018..c140d98ae 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 2c7b7ecfb..c08729022 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 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: diff --git a/src/surface.cpp b/src/surface.cpp index d60b0c52e..9323629cd 100644 --- a/src/surface.cpp +++ b/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;