From 00b96bba7841ed5da179a1b66c30f55355c504f2 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Fri, 27 Jul 2012 14:38:24 -0400 Subject: [PATCH] Added description of neighbor lists and derivations for reflective boundary conditions. --- docs/source/methods/geometry.rst | 193 +++++++++++++++++++++++++++++-- 1 file changed, 183 insertions(+), 10 deletions(-) diff --git a/docs/source/methods/geometry.rst b/docs/source/methods/geometry.rst index 7dabd81be3..b13efbfe20 100644 --- a/docs/source/methods/geometry.rst +++ b/docs/source/methods/geometry.rst @@ -357,26 +357,199 @@ search is done over every cell in the base universe. Building Neighbor Lists ----------------------- +After the geometry has been loaded and stored in memory from an input file, +OpenMC builts a list for each surface containing any cells that contain the +surface in their specification in order to speed up processing of surface +crossings. The algorithm to build these lists is as follows. First, we loop over +all cells in the geometry and count up how many times each surface appears in a +specification as bounding a negative half-space and bounding a positive +half-space. Two arrays are then allocated for each surface, one that lists each +cell that contains the negative half-space of the surface and one that lists +each cell that contains the positive half-space of the surface. Another loop is +performed over all cells and the neighbor lists are populated for each surface. + .. _reflection: -------------------- -Reflective Surfaces -------------------- +------------------------------ +Reflective Boundary Conditions +------------------------------ -In general, a surface can be written in the form :math:`f(x,y,z) = 0`. If a -neutron is traveling in direction :math:`\mathbf{v}` and crosses a reflective -surface of the above form, it can be shown that the velocity vector will then -become +If the velocity of a particle is :math:`\mathbf{v}` and it crosses a surface of +the form :math:`f(x,y,z) = 0` with a reflective boundary condition, it can be +shown based on geometric arguments that the velocity vector will then become .. math:: + :label: reflection-v \mathbf{v'} = \mathbf{v} - 2 (\mathbf{v} \cdot \hat{\mathbf{n}}) \hat{\mathbf{n}} where :math:`\hat{\mathbf{n}}` is a unit vector normal to the surface at the -point of the surface crossing. The direction of the surface normal will be the -gradient to the surface at the point of crossing, i.e. :math:`\mathbf{n} = -\nabla f(x,y,z)`. +point of the surface crossing. The rationale for this can be understood by +noting that :math:`(\mathbf{v} \cdot \hat{\mathbf{n}}) \hat{\mathbf{n}}` is the +projection of the velocity vector onto the normal vector. By subtracting two +times this projection, the velocity is reflected with respect to the surface +normal. Since the velocity of the particle will not change as it undergoes +reflection, we can work with the direction of the particle instead, simplifying +equation :eq:`reflection-v` to + +.. math:: + :label: reflection-omega + + \mathbf{\Omega'} = \mathbf{\Omega} - 2 (\mathbf{\Omega} \cdot + \hat{\mathbf{n}}) \hat{\mathbf{n}} + + +The direction of the surface normal will be the gradient to the surface at the +point of crossing, i.e. :math:`\mathbf{n} = \nabla f(x,y,z)`. Substituting this +into equation :eq:`reflection-omega`, we get + +.. math:: + :label: reflection-omega-2 + + \mathbf{\Omega'} = \mathbf{\Omega} - \frac{2 ( \mathbf{\Omega} \cdot \nabla + f )}{|| \nabla f ||^2} \nabla f + + +If we write the initial and final directions in terms of their vector +components, :math:`\mathbf{\Omega} = (u,v,w)` and :math:`\mathbf{\Omega'} = (u', +v', w')`, this allows us to represent equation :eq:`reflection-omega` as a +series of equations: + +.. math:: + :label: reflection-system + + u' = u - \frac{2 ( \mathbf{\Omega} \cdot \nabla f )}{|| \nabla f ||^2} + \frac{\partial f}{\partial x} \\ + + v' = v - \frac{2 ( \mathbf{\Omega} \cdot \nabla f )}{|| \nabla f ||^2} + \frac{\partial f}{\partial y} \\ + + w' = w - \frac{2 ( \mathbf{\Omega} \cdot \nabla f )}{|| \nabla f ||^2} + \frac{\partial f}{\partial z} + +We can now use this form to develop rules for how to transform a particle's +direction for different types of surfaces. + +Plane Perpendicular to an Axis +------------------------------ + +For a plane that is perpendicular to an axis, the rule for reflection is almost +so simple that no derivation is needed at all. Nevertheless, we will proceed +with the derivation to confirm that the rules of geometry agree with our +intuition. The gradient of the surface :math:`f(x,y,z) = x - x_0 = 0` is simply +:math:`\nabla f = (1, 0, 0)`. Note that this vector is already normalized, +i.e. :math:`|| \nabla f || = 1`. The second two equations in +:eq:`reflection-system` tell us that :math:`v` and :math:`w` do not change and +the first tell us that + +.. math:: + :label: reflection-xplane + + u' = u - 2u = -u + +We see that reflection for a plane perpendicular to an axis only entails +negating the directional cosine for that axis. + +Generic Plane +------------- + +A generic plane has the form :math:`f(x,y,z) = Ax + By + Cz - D = 0`. Thus, the +gradient to the surface is simply :math:`\nabla f = (A,B,C)` whose norm squared +is :math:`A^2 + B^2 + C^2`. This implies that + +.. math:: + :label: reflection-plane-constant + + \frac{2 (\mathbf{\Omega} \cdot \nabla f)}{|| \nabla f ||^2} = \frac{2(Au + + Bv + Cw)}{A^2 + B^2 + C^2} + +Substituting equation :eq:`reflection-plane-constant` into equation +:eq:`reflection-system` gives us the form of the solution. For example, the +x-component of the reflected direction will be + +.. math:: + :label: reflection-plane + + u' = u - \frac{2A(Au + Bv + Cw)}{A^2 + B^2 + C^2} + + +Cylinder Parallel to an Axis +---------------------------- + +A cylinder parallel to, for example, the x-axis has the form :math:`f(x,y,z) = +(y - y_0)^2 + (z - z_0)^2 - R^2 = 0`. Thus, the gradient to the surface is + +.. math:: :label: reflection-cylinder-grad + + \nabla f = 2 \left ( \begin{array}{c} 0 \\ y - y_0 \\ z - z_0 \end{array} + \right ) = 2 \left ( \begin{array}{c} 0 \\ \bar{y} \\ \bar{z} \end{array} + \right ) + +where we have introduced the constants :math:`\bar{y}` and +:math:`\bar{z}`. Taking the square of the norm of the gradient, we find that + +.. math:: :label: reflection-cylinder-norm + + || \nabla f ||^2 = 4 \bar{y}^2 + 4 \bar{z}^2 = 4 R^2 + +This imples that + +.. math:: :label: reflection-cylinder-constant + + \frac{2 (\mathbf{\Omega} \cdot \nabla f)}{|| \nabla f ||^2} = + \frac{\bar{y}v + \bar{z}w}{R^2} + +Substituting equations :eq:`reflection-cylinder-constant` and +:eq:`reflection-cylinder-grad` into equation :eq:`reflection-system` gives us +the form of the solution. In this case, the x-component will not change. The y- +and z-components of the reflected direction will be + +.. math:: :label: reflection-cylinder + + v' = v - \frac{2 ( \bar{y}v + \bar{z}w ) \bar{y}}{R^2} \\ + + w' = w - \frac{2 ( \bar{y}v + \bar{z}w ) \bar{z}}{R^2} + + +Sphere +------ + +The surface equation for a sphere has the form :math:`f(x,y,z) = (x - x_0)^2 + +(y - y_0)^2 + (z - z_0)^2 - R^2 = 0`. Thus, the gradient to the surface is + +.. math:: :label: reflection-sphere-grad + + \nabla f = 2 \left ( \begin{array}{c} x - x_0 \\ y - y_0 \\ z - z_0 + \end{array} \right ) = 2 \left ( \begin{array}{c} \bar{x} \\ \bar{y} \\ + \bar{z} \end{array} \right ) + +where we have introduced the constants :math:`\bar{x}, \bar{y}, \bar{z}`. Taking +the square of the norm of the gradient, we find that + +.. math:: :label: reflection-sphere-norm + + || \nabla f ||^2 = 4 \bar{x}^2 + 4 \bar{y}^2 + 4 \bar{z}^2 = 4 R^2 + +This imples that + +.. math:: :label: reflection-sphere-constant + + \frac{2 (\mathbf{\Omega} \cdot \nabla f)}{|| \nabla f ||^2} = + \frac{\bar{x}u + \bar{y}v + \bar{z}w}{R^2} + +Substituting equations :eq:`reflection-sphere-constant` and +:eq:`reflection-sphere-grad` into equation :eq:`reflection-system` gives us the +form of the solution: + +.. math:: :label: reflection-sphere + + u' = u - \frac{2 ( \bar{x}u + \bar{y}v + \bar{z}w ) \bar{x} }{R^2} \\ + + v' = v - \frac{2 ( \bar{x}u + \bar{y}v + \bar{z}w ) \bar{y} }{R^2} \\ + + w' = w - \frac{2 ( \bar{x}u + \bar{y}v + \bar{z}w ) \bar{z} }{R^2} \\ + .. _constructive solid geometry: http://en.wikipedia.org/wiki/Constructive_solid_geometry .. _surfaces: http://en.wikipedia.org/wiki/Surface