Merge remote-tracking branch 'upstream/develop' into velocity-tally

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/statepoint.py

@@ -358,7 +358,7 @@ class StatePoint(object):
 
                 # Create Tally object and assign basic properties
                 tally = openmc.Tally(tally_id=tally_key)
-                tally._statepoint = self
+                tally._sp_filename = self._f.filename
                 tally.estimator = self._f['{0}{1}/estimator'.format(
                     base, tally_key)].value.decode()
                 tally.num_realizations = n_realizations

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.
 

openmc/tallies.py

@@ -103,7 +103,7 @@ class Tally(object):
         self._with_batch_statistics = False
         self._derived = False
 
-        self._statepoint = None
+        self._sp_filename = None
         self._results_read = False
 
     def __deepcopy__(self, memo):
@@ -124,7 +124,7 @@ class Tally(object):
             clone._with_summary = self.with_summary
             clone._with_batch_statistics = self.with_batch_statistics
             clone._derived = self.derived
-            clone._statepoint = self._statepoint
+            clone._sp_filename = self._sp_filename
             clone._results_read = self._results_read
 
             clone._filters = []
@@ -264,12 +264,17 @@ class Tally(object):
 
     @property
     def sum(self):
-        if not self._statepoint:
+        if not self._sp_filename:
             return None
 
         if not self._results_read:
+            import h5py
+
+            # Open the HDF5 statepoint file
+            f = h5py.File(self._sp_filename, 'r')
+
             # Extract Tally data from the file
-            data = self._statepoint._f['tallies/tally {0}/results'.format(
+            data = f['tallies/tally {0}/results'.format(
                 self.id)].value
             sum = data['sum']
             sum_sq = data['sum_sq']
@@ -293,11 +298,14 @@ class Tally(object):
             # Indicate that Tally results have been read
             self._results_read = True
 
+            # Close the HDF5 statepoint file
+            f.close()
+
         return self._sum
 
     @property
     def sum_sq(self):
-        if not self._statepoint:
+        if not self._sp_filename:
             return None
 
         if not self._results_read:
@@ -309,7 +317,7 @@ class Tally(object):
     @property
     def mean(self):
         if self._mean is None:
-            if not self._statepoint:
+            if not self._sp_filename:
                 return None
 
             self._mean = self.sum / self.num_realizations
@@ -318,7 +326,7 @@ class Tally(object):
     @property
     def std_dev(self):
         if self._std_dev is None:
-            if not self._statepoint:
+            if not self._sp_filename:
                 return None
 
             n = self.num_realizations
@@ -426,8 +434,13 @@ class Tally(object):
         # If the score is already in the Tally, don't add it again
         if score in self.scores:
             return
-        else:
+
+        # Normal score strings
+        if isinstance(score, basestring):
             self._scores.append(score.strip())
+        # CrossScores
+        else:
+            self._scores.append(score)
 
     @num_score_bins.setter
     def num_score_bins(self, num_score_bins):

src/ace.F90

@@ -215,6 +215,16 @@ contains
 
     end do MATERIAL_LOOP3
 
+    ! Show which nuclide results in lowest energy for neutron transport
+    do i = 1, n_nuclides_total
+      if (nuclides(i)%energy(nuclides(i)%n_grid) == energy_max_neutron) then
+        call write_message("Maximum neutron transport energy: " // &
+             trim(to_str(energy_max_neutron)) // " MeV for " // &
+             trim(adjustl(nuclides(i)%name)), 6)
+        exit
+      end if
+    end do
+
   end subroutine read_xs
 
 !===============================================================================
@@ -482,6 +492,10 @@ contains
       ! Continue reading elastic scattering and heating
       nuc % elastic = get_real(NE)
 
+      ! Determine if minimum/maximum energy for this nuclide is greater/less
+      ! than the previous
+      energy_min_neutron = max(energy_min_neutron, nuc%energy(1))
+      energy_max_neutron = min(energy_max_neutron, nuc%energy(NE))
     end if
 
   end subroutine read_esz

src/cross_section.F90

@@ -56,7 +56,7 @@ contains
     if (grid_method == GRID_MAT_UNION) then
       call find_energy_index(p % E, p % material)
     else if (grid_method == GRID_LOGARITHM) then
-      u = int(log(p % E/1.0e-11_8)/log_spacing)
+      u = int(log(p % E/energy_min_neutron)/log_spacing)
     end if
 
     ! Determine if this material has S(a,b) tables

src/energy_grid.F90

@@ -73,8 +73,8 @@ contains
     type(Nuclide), pointer :: nuc
 
     ! Set minimum/maximum energies
-    E_max = 20.0_8
-    E_min = 1.0e-11_8
+    E_max = energy_max_neutron
+    E_min = energy_min_neutron
 
     ! Determine equal-logarithmic energy spacing
     M = n_log_bins

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/global.F90

@@ -74,6 +74,10 @@ module global
   integer :: n_sab_tables     ! Number of S(a,b) thermal scattering tables
   integer :: n_listings       ! Number of listings in cross_sections.xml
 
+  ! Minimum/maximum energies
+  real(8) :: energy_min_neutron = ZERO
+  real(8) :: energy_max_neutron = INFINITY
+
   ! Dictionaries to look up cross sections and listings
   type(DictCharInt) :: nuclide_dict
   type(DictCharInt) :: sab_dict

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
@@ -2847,6 +2861,10 @@ contains
           ! MOMENT_STRS(:)
           ! If so, check the order, store if OK, then reset the number to 'n'
           score_name = trim(sarray(j))
+
+          ! Append the score to the list of possible trigger scores
+          if (trigger_on) call trigger_scores % add_key(trim(score_name), j)
+
           do imomstr = 1, size(MOMENT_STRS)
             if (starts_with(score_name,trim(MOMENT_STRS(imomstr)))) then
               n_order_pos = scan(score_name,'0123456789')
@@ -3263,11 +3281,8 @@ contains
             end if
 
           end select
-
-          ! Append the score to the list of possible trigger scores
-          if (trigger_on) call trigger_scores % add_key(trim(score_name), l)
-
         end do
+
         t % n_score_bins = n_scores
         t % n_user_score_bins = n_words
 

src/physics.F90

@@ -1250,8 +1250,8 @@ contains
           call sample_energy(edist, p % E, E_out)
         end if
 
-        ! resample if energy is >= 20 MeV
-        if (E_out < 20) exit
+        ! resample if energy is greater than maximum neutron energy
+        if (E_out < energy_max_neutron) exit
 
         ! check for large number of resamples
         n_sample = n_sample + 1
@@ -1279,8 +1279,8 @@ contains
           call sample_energy(rxn%edist, p % E, E_out)
         end if
 
-        ! resample if energy is >= 20 MeV
-        if (E_out < 20) exit
+        ! resample if energy is greater than maximum neutron energy
+        if (E_out < energy_max_neutron) exit
 
         ! check for large number of resamples
         n_sample = n_sample + 1

src/source.F90

@@ -211,8 +211,9 @@ contains
     case (SRC_ENERGY_MONO)
       ! Monoenergtic source
       site%E = external_source%params_energy(1)
-      if (site%E >= 20) then
-        call fatal_error("Source energies above 20 MeV not allowed.")
+      if (site%E >= energy_max_neutron) then
+        call fatal_error("Source energy above range of energies of at least &
+             &one cross section table")
       end if
 
     case (SRC_ENERGY_MAXWELL)
@@ -221,8 +222,8 @@ contains
         ! Sample Maxwellian fission spectrum
         site%E = maxwell_spectrum(a)
 
-        ! resample if energy is >= 20 MeV
-        if (site%E < 20) exit
+        ! resample if energy is greater than maximum neutron energy
+        if (site%E < energy_max_neutron) exit
       end do
 
     case (SRC_ENERGY_WATT)
@@ -232,8 +233,8 @@ contains
         ! Sample Watt fission spectrum
         site%E = watt_spectrum(a, b)
 
-        ! resample if energy is >= 20 MeV
-        if (site%E < 20) exit
+        ! resample if energy is greater than maximum neutron energy
+        if (site%E < energy_max_neutron) exit
       end do
 
     case default

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

src/trigger.F90

@@ -92,7 +92,6 @@ contains
     character(len=52), intent(inout) :: name  ! "eigenvalue" or tally score
 
     integer :: i              ! index in tallies array
-    integer :: j              ! level in tally hierarchy
     integer :: n              ! loop index for nuclides
     integer :: s              ! loop index for triggers
     integer :: filter_index   ! index in results array for filters
@@ -167,28 +166,8 @@ contains
 
             ! Initialize bins, filter level
             matching_bins(1:t % n_filters) = 0
-            j = 1
 
-            ! Find filter index
-            FILTER_LOOP: do
-              find_bin: do
-                if (t % n_filters == 0) exit find_bin
-                matching_bins(j) = matching_bins(j) + 1
-                if (matching_bins(j) > t % filters(j) % n_bins) then
-                  if (j == 1) exit FILTER_LOOP
-                  matching_bins(j) = 0
-                  j = j - 1
-                else
-                  if (j == t % n_filters) exit find_bin
-                end if
-              end do find_bin
-
-              if (t % n_filters > 0) then
-                filter_index = sum((max(matching_bins(1:t%n_filters),1) - 1) * &
-                     t % stride) + 1
-              else
-                filter_index = 1
-              end if
+            FILTER_LOOP: do filter_index = 1, t % total_filter_bins
 
               ! Initialize score index
               score_index = trigger % score_index

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()