Merge pull request #477 from paulromano/quadric

General quadric surface support

docs/source/methods/geometry.rst

@@ -10,7 +10,7 @@ Constructive Solid Geometry
 
 OpenMC uses a technique known as `constructive solid geometry`_ (CSG) to build
 arbitrarily complex three-dimensional models in Euclidean space. In a CSG model,
-every unique object is described as the union, intersection, or difference of
+every unique object is described as the union and/or intersection of
 *half-spaces* created by bounding `surfaces`_. Every surface divides all of
 space into exactly two half-spaces. We can mathematically define a surface as a
 collection of points that satisfy an equation of the form :math:`f(x,y,z) = 0`
@@ -54,13 +54,12 @@ dividing space into two half-spaces.
    Example of an ellipse and its associated half-spaces.
 
 References to half-spaces created by surfaces are used to define regions of
-space of uniform composition, known as *cells*. While some codes allow regions
-to be defined by intersections, unions, and differences or half-spaces, OpenMC
-is currently limited to cells defined only as intersections of
-half-spaces. Thus, the specification of the cell must include a list of
-half-space references whose intersection defines the region. The region is then
-assigned a material defined elsewhere. Figure :num:`fig-union` shows an
-example of a cell defined as the intersection of an ellipse and two planes.
+space of uniform composition, which are then assigned to *cells*. OpenMC allows
+regions to be defined using union, intersection, and complement operators. As in
+MCNP_, the intersection operator is implicit as doesn't need to be written in a
+region specification. A defined region is then associated with a material
+composition in a cell. Figure :num:`fig-union` shows an example of a cell region
+defined as the intersection of an ellipse and two planes.
 
 .. _fig-union:
 
@@ -117,6 +116,10 @@ to fully define the surface.
     | Cone parallel to the | z-cone     | :math:`(x-x_0)^2 + (y-y_0)^2 | :math:`x_0 \; y_0 \;    |
     | :math:`z`-axis       |            | = R^2(z-z_0)^2`              | z_0 \; R^2`             |
     +----------------------+------------+------------------------------+-------------------------+
+    | 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`                 |
+    +----------------------+------------+------------------------------+-------------------------+
 
 .. _universes:
 

docs/source/usersguide/input.rst

@@ -787,7 +787,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", or "z-cone".
+    "y-cone", "z-cone", or "quadric".
 
     *Default*: None
 
@@ -855,6 +855,12 @@ The following quadratic surfaces can be modeled:
     R^2 (z - z_0)^2`. The coefficients specified are ":math:`x_0 \: y_0 \: z_0
     \: R^2`".
 
+  :quadric:
+     A general quadric surface of the form :math:`Ax^2 + By^2 + Cz^2 + Dxy +
+     Eyz + Fxz + Gx + Hy + Jz + K = 0` The coefficients specified are ":math:`A
+     \: B \: C \: D \: E \: F \: G \: H \: J \: K`".
+
+
 ``<cell>`` Element
 ------------------
 

docs/source/usersguide/output/summary.rst

@@ -132,7 +132,8 @@ The current revision of the summary file format is 1.
 **/geometry/surfaces/surface <uid>/type** (*char[]*)
 
     Type of the surface. Can be 'x-plane', 'y-plane', 'z-plane', 'plane',
-    'x-cylinder', 'y-cylinder', 'sphere', 'x-cone', 'y-cone', or 'z-cone'.
+    'x-cylinder', 'y-cylinder', 'sphere', 'x-cone', 'y-cone', 'z-cone', or
+    'quadric'.
 
 **/geometry/surfaces/surface <uid>/coefficients** (*double[]*)
 

openmc/summary.py

@@ -179,6 +179,11 @@ class Summary(object):
                 if surf_type == 'z-cone':
                     surface = openmc.ZCone(surface_id, bc, x0, y0, z0, R2, name)
 
+            elif surf_type == 'quadric':
+                a, b, c, d, e, f, g, h, j, k = coeffs
+                surface = openmc.Quadric(surface_id, bc, a, b, c, d, e, f,
+                                         g, h, j, k, name)
+
             # Add Surface to global dictionary of all Surfaces
             self.surfaces[index] = surface
 

openmc/surface.py

@@ -949,6 +949,152 @@ class ZCone(Cone):
         self._type = 'z-cone'
 
 
+class Quadric(Surface):
+    """A sphere of the form :math:`Ax^2 + By^2 + Cz^2 + Dxy + Eyz + Fxz + Gx + Hy +
+    Jz + K`.
+
+    Parameters
+    ----------
+    surface_id : int
+        Unique identifier for the surface. If not specified, an identifier will
+        automatically be assigned.
+    boundary_type : {'transmission, 'vacuum', 'reflective', 'periodic'}
+        Boundary condition that defines the behavior for particles hitting the
+        surface. Defaults to transmissive boundary condition where particles
+        freely pass through the surface.
+    a, b, c, d, e, f, g, h, j, k : float
+        coefficients for the surface
+    name : str
+        Name of the sphere. If not specified, the name will be the empty string.
+
+    Attributes
+    ----------
+    a, b, c, d, e, f, g, h, j, k : float
+        coefficients for the surface
+
+    """
+
+    def __init__(self, surface_id=None, boundary_type='transmission',
+                 a=None, b=None, c=None, d=None, e=None, f=None, g=None,
+                 h=None, j=None, k=None, name=''):
+        # Initialize Quadric class attributes
+        super(Quadric, self).__init__(surface_id, boundary_type, name=name)
+
+        self._type = 'quadric'
+        self._coeff_keys = ['a', 'b', 'c', 'd', 'e', 'f', 'g', 'h', 'j', 'k']
+
+        if a is not None:
+            self.a = a
+        if b is not None:
+            self.b = b
+        if c is not None:
+            self.c = c
+        if d is not None:
+            self.d = d
+        if e is not None:
+            self.e = e
+        if f is not None:
+            self.f = f
+        if g is not None:
+            self.g = g
+        if h is not None:
+            self.h = h
+        if j is not None:
+            self.j = j
+        if k is not None:
+            self.k = k
+
+    @property
+    def a(self):
+        return self.coeffs['a']
+
+    @property
+    def b(self):
+        return self.coeffs['b']
+
+    @property
+    def c(self):
+        return self.coeffs['c']
+
+    @property
+    def d(self):
+        return self.coeffs['d']
+
+    @property
+    def e(self):
+        return self.coeffs['e']
+
+    @property
+    def f(self):
+        return self.coeffs['f']
+
+    @property
+    def g(self):
+        return self.coeffs['g']
+
+    @property
+    def h(self):
+        return self.coeffs['h']
+
+    @property
+    def j(self):
+        return self.coeffs['j']
+
+    @property
+    def k(self):
+        return self.coeffs['k']
+
+    @a.setter
+    def a(self, a):
+        check_type('a coefficient', a, Real)
+        self._coeffs['a'] = a
+
+    @b.setter
+    def b(self, b):
+        check_type('b coefficient', b, Real)
+        self._coeffs['b'] = b
+
+    @c.setter
+    def c(self, c):
+        check_type('c coefficient', c, Real)
+        self._coeffs['c'] = c
+
+    @d.setter
+    def d(self, d):
+        check_type('d coefficient', d, Real)
+        self._coeffs['d'] = d
+
+    @e.setter
+    def e(self, e):
+        check_type('e coefficient', e, Real)
+        self._coeffs['e'] = e
+
+    @f.setter
+    def f(self, f):
+        check_type('f coefficient', f, Real)
+        self._coeffs['f'] = f
+
+    @g.setter
+    def g(self, g):
+        check_type('g coefficient', g, Real)
+        self._coeffs['g'] = g
+
+    @h.setter
+    def h(self, h):
+        check_type('h coefficient', h, Real)
+        self._coeffs['h'] = h
+
+    @j.setter
+    def j(self, j):
+        check_type('j coefficient', j, Real)
+        self._coeffs['j'] = j
+
+    @k.setter
+    def k(self, k):
+        check_type('k coefficient', k, Real)
+        self._coeffs['k'] = k
+
+
 class Halfspace(Region):
     """A positive or negative half-space region.
 

src/geometry.F90

@@ -430,7 +430,7 @@ contains
       call find_cell(p, found)
       if (.not. found) then
         call handle_lost_particle(p, "Couldn't find particle after reflecting&
-             & from surface.")
+             & from surface " // trim(to_str(surf%id)) // ".")
         return
       end if
 

src/input_xml.F90

@@ -1294,6 +1294,9 @@ contains
       case ('z-cone')
         coeffs_reqd  = 4
         allocate(SurfaceZCone :: surfaces(i)%obj)
+      case ('quadric')
+        coeffs_reqd  = 10
+        allocate(SurfaceQuadric :: surfaces(i)%obj)
       case default
         call fatal_error("Invalid surface type: " // trim(word))
       end select
@@ -1378,6 +1381,17 @@ contains
         s%y0 = coeffs(2)
         s%z0 = coeffs(3)
         s%r2 = coeffs(4)
+      type is (SurfaceQuadric)
+        s%A = coeffs(1)
+        s%B = coeffs(2)
+        s%C = coeffs(3)
+        s%D = coeffs(4)
+        s%E = coeffs(5)
+        s%F = coeffs(6)
+        s%G = coeffs(7)
+        s%H = coeffs(8)
+        s%J = coeffs(9)
+        s%K = coeffs(10)
       end select
 
       ! No longer need coefficients

src/summary.F90

@@ -276,6 +276,11 @@ contains
         allocate(coeffs(4))
         coeffs(:) = [s%x0, s%y0, s%z0, s%r2]
 
+      type is (SurfaceQuadric)
+        call write_dataset(surface_group, "type", "quadric")
+        allocate(coeffs(10))
+        coeffs(:) = [s%A, s%B, s%C, s%D, s%E, s%F, s%G, s%H, s%J, s%K]
+
       end select
       call write_dataset(surface_group, "coefficients", coeffs)
       deallocate(coeffs)

src/surface_header.F90

@@ -1,6 +1,6 @@
 module surface_header
 
-  use constants, only: ONE, TWO, ZERO, INFINITY, FP_COINCIDENT
+  use constants, only: ONE, TWO, ZERO, HALF, INFINITY, FP_COINCIDENT
 
   implicit none
 
@@ -183,6 +183,15 @@ module surface_header
     procedure :: normal => z_cone_normal
   end type SurfaceZCone
 
+  type, extends(Surface) :: SurfaceQuadric
+    ! Ax^2 + By^2 + Cz^2 + Dxy + Eyz + Fxz + Gx + Hy + Jz + K = 0
+    real(8) :: A, B, C, D, E, F, G, H, J, K
+  contains
+    procedure :: evaluate => quadric_evaluate
+    procedure :: distance => quadric_distance
+    procedure :: normal => quadric_normal
+  end type SurfaceQuadric
+
 contains
 
 !===============================================================================
@@ -640,7 +649,7 @@ contains
   end function z_cylinder_normal
 
 !===============================================================================
-! SphereImplementation
+! SurfaceSphere Implementation
 !===============================================================================
 
   pure function sphere_evaluate(this, xyz) result(f)
@@ -874,7 +883,7 @@ contains
   end function y_cone_normal
 
 !===============================================================================
-! SurfaceZConeImplementation
+! SurfaceZCone Implementation
 !===============================================================================
 
   pure function z_cone_evaluate(this, xyz) result(f)
@@ -953,4 +962,91 @@ contains
     uvw(3) = -TWO*this%r2*(xyz(3) - this%z0)
   end function z_cone_normal
 
+!===============================================================================
+! SurfaceQuadric Implementation
+!===============================================================================
+
+  pure function quadric_evaluate(this, xyz) result(f)
+    class(SurfaceQuadric), intent(in) :: this
+    real(8), intent(in) :: xyz(3)
+    real(8) :: f
+
+    associate (x => xyz(1), y => xyz(2), z => xyz(3))
+      f = x*(this%A*x + this%D*y + this%G) + &
+           y*(this%B*y + this%E*z + this%H) + &
+           z*(this%C*z + this%F*x + this%J) + this%K
+    end associate
+  end function quadric_evaluate
+
+  pure function quadric_distance(this, xyz, uvw, coincident) result(d)
+    class(SurfaceQuadric), intent(in) :: this
+    real(8), intent(in) :: xyz(3)
+    real(8), intent(in) :: uvw(3)
+    logical, intent(in) :: coincident
+    real(8) :: d
+
+    real(8) :: a, k, c
+    real(8) :: quad, b
+
+    associate (x => xyz(1), y => xyz(2), z => xyz(3), &
+         u => uvw(1), v => uvw(2), w => uvw(3))
+
+      a = this%A*u*u + this%B*v*v + this%C*w*w + this%D*u*v + this%E*v*w + &
+           this%F*u*w
+      k = (this%A*u*x + this%B*v*y + this%C*w*z + HALF*(this%D*(u*y + v*x) + &
+           this%E*(v*z + w*y) + this%F*(w*x + u*z) + this%G*u + this%H*v + &
+           this%J*w))
+      c = this%A*x*x + this%B*y*y + this%C*z*z + this%D*x*y + this%E*y*z + &
+           this%F*x*z + this%G*x + this%H*y + this%J*z + this%K
+      quad = k*k - a*c
+
+      if (quad < ZERO) then
+        ! no intersection with surface
+
+        d = INFINITY
+
+      elseif (coincident .or. abs(c) < FP_COINCIDENT) then
+        ! particle is on the surface, thus one distance is positive/negative and
+        ! the other is zero. The sign of k determines which distance is zero and
+        ! which is not.
+
+        if (k >= ZERO) then
+          d = (-k - sqrt(quad))/a
+        else
+          d = (-k + sqrt(quad))/a
+        end if
+
+      else
+        ! calculate both solutions to the quadratic
+        quad = sqrt(quad)
+        d = (-k - quad)/a
+        b = (-k + quad)/a
+
+        ! determine the smallest positive solution
+        if (d < ZERO) then
+          if (b > ZERO) then
+            d = b
+          end if
+        else
+          if (b > ZERO) d = min(d, b)
+        end if
+      end if
+
+      ! If the distance was negative, set boundary distance to infinity
+      if (d <= ZERO) d = INFINITY
+    end associate
+  end function quadric_distance
+
+  pure function quadric_normal(this, xyz) result(uvw)
+    class(SurfaceQuadric), intent(in) :: this
+    real(8), intent(in) :: xyz(3)
+    real(8) :: uvw(3)
+
+    associate (x => xyz(1), y => xyz(2), z => xyz(3))
+      uvw(1) = TWO*this%A*x + this%D*y + this%F*z + this%G
+      uvw(2) = TWO*this%B*y + this%D*x + this%E*z + this%H
+      uvw(3) = TWO*this%C*z + this%E*y + this%F*x + this%J
+    end associate
+  end function quadric_normal
+
 end module surface_header

tests/test_quadric_surfaces/geometry.xml

@@ -0,0 +1,14 @@
+<?xml version="1.0"?>
+<geometry>
+
+  <surface id="1" type="sphere" coeffs="0. 0. 5. 5." boundary="reflective" />
+  <surface id="2" type="quadric" coeffs="1. 1. 1. 0. 0. 0. 0. 0. 0. -81." boundary="reflective" />
+  <surface id="3" type="z-plane" coeffs="5." />
+  <surface id="4" type="z-cylinder" coeffs="0. 0. 5." boundary="reflective" />
+  <surface id="5" type="z-cone" coeffs="0. 0. -10. 1." boundary="reflective" />
+  <surface id="6" type="plane" coeffs="0.2 0.2 1.0 -8." boundary="reflective" />
+
+  <cell id="1" material="1" region="-1 -2 3" />
+  <cell id="2" material="1" region="-3 -4 -5 6" />
+
+</geometry>

tests/test_quadric_surfaces/materials.xml

@@ -3,7 +3,8 @@
 
   <material id="1">
     <density value="4.5" units="g/cc" />
-    <nuclide name="U-235" xs="71c" ao="1.0" />
+    <nuclide name="U-238" xs="71c" ao="1.0" />
+    <nuclide name="U-235" xs="71c" ao="0.06" />
   </material>
 
 </materials>

tests/test_quadric_surfaces/results_true.dat

@@ -0,0 +1,2 @@
+k-combined:
+9.706301E-01 4.351374E-02

tests/test_quadric_surfaces/settings.xml

@@ -8,7 +8,7 @@
   </eigenvalue>
 
   <source>
-    <space type="point" parameters="0 0 -5" />
+    <space type="point" parameters="0. 0. 0." />
   </source>
 
 </settings>

tests/test_reflective_cone/geometry.xml

@@ -1,7 +0,0 @@
-<?xml version="1.0"?>
-<geometry>
-
-  <surface id="1" type="z-cone" coeffs="0 0 0 5" boundary="reflective"/>
-  <cell id="1" material="1" region="-1" />
-    
-</geometry>

tests/test_reflective_cone/materials.xml

@@ -1,9 +0,0 @@
-<?xml version="1.0"?>
-<materials>
-
-  <material id="1">
-    <density value="4.5" units="g/cc" />
-    <nuclide name="U-235" xs="71c" ao="1.0" />
-  </material>
-
-</materials>

tests/test_reflective_cone/results_true.dat

@@ -1,2 +0,0 @@
-k-combined:
-2.269987E+00 4.469683E-03

tests/test_reflective_cylinder/geometry.xml

@@ -1,8 +0,0 @@
-<?xml version="1.0"?>
-<geometry>
-
-  <!-- cylinder with radius 10 -->
-  <surface id="1" type="z-cylinder" coeffs="0 0 10" boundary="reflective"/>
-  <cell id="1" material="1" region="-1" />
-    
-</geometry>

tests/test_reflective_cylinder/materials.xml

@@ -1,9 +0,0 @@
-<?xml version="1.0"?>
-<materials>
-
-  <material id="1">
-    <density value="4.5" units="g/cc" />
-    <nuclide name="U-235" xs="71c" ao="1.0" />
-  </material>
-
-</materials>

tests/test_reflective_cylinder/results_true.dat

@@ -1,2 +0,0 @@
-k-combined:
-2.272436E+00 7.831006E-04

tests/test_reflective_cylinder/settings.xml

@@ -1,16 +0,0 @@
-<?xml version="1.0"?>
-<settings>
-
-  <eigenvalue>
-    <batches>10</batches>
-    <inactive>5</inactive>
-    <particles>1000</particles>
-  </eigenvalue>
-
-  <source>
-    <space type="box">
-      <parameters>-4 -4 -4  4  4  4</parameters>
-    </space>
-  </source>
-
-</settings>

tests/test_reflective_cylinder/test_reflective_cylinder.py

@@ -1,11 +0,0 @@
-#!/usr/bin/env python
-
-import os
-import sys
-sys.path.insert(0, os.pardir)
-from testing_harness import TestHarness
-
-
-if __name__ == '__main__':
-    harness = TestHarness('statepoint.10.*')
-    harness.main()

tests/test_reflective_sphere/geometry.xml

@@ -1,8 +0,0 @@
-<?xml version="1.0"?>
-<geometry>
-
-  <!-- Sphere with radius 10 -->
-  <surface id="1" type="sphere" coeffs="0 0 0 10" boundary="reflective"/>
-  <cell id="1" material="1" region="-1" />
-    
-</geometry>

tests/test_reflective_sphere/results_true.dat

@@ -1,2 +0,0 @@
-k-combined:
-2.271012E+00 3.466351E-03

tests/test_reflective_sphere/settings.xml

@@ -1,16 +0,0 @@
-<?xml version="1.0"?>
-<settings>
-
-  <eigenvalue>
-    <batches>10</batches>
-    <inactive>5</inactive>
-    <particles>1000</particles>
-  </eigenvalue>
-
-  <source>
-    <space type="box">
-      <parameters>-4 -4 -4  4  4  4</parameters>
-    </space>
-  </source>
-
-</settings>

tests/test_reflective_sphere/test_reflective_sphere.py

@@ -1,11 +0,0 @@
-#!/usr/bin/env python
-
-import os
-import sys
-sys.path.insert(0, os.pardir)
-from testing_harness import TestHarness
-
-
-if __name__ == '__main__':
-    harness = TestHarness('statepoint.10.*')
-    harness.main()