Merge pull request #1289 from paulromano/cylinder-fix-again

Fix cylinder_from_points (again)

docs/source/io_formats/source.rst

@@ -13,8 +13,9 @@ is that documented here.
 :Attributes: - **filetype** (*char[]*) -- String indicating the type of file.
 
 :Datasets:
+
            - **source_bank** (Compound type) -- Source bank information for each
              particle. The compound type has fields ``wgt``, ``xyz``, ``uvw``,
-             ``E``, and ``delayed_group``, which represent the weight, position,
-             direction, energy, energy group, and delayed_group of the source
-             particle, respectively.
+             ``E``, ``delayed_group``, and ``particle``, which represent the
+             weight, position, direction, energy, energy group, delayed group,
+             and type of the source particle, respectively.

openmc/data/neutron.py

@@ -801,15 +801,14 @@ class IncidentNeutron(EqualityMixin):
         """
         with tempfile.TemporaryDirectory() as tmpdir:
             # Run NJOY to create an ACE library
-            ace_file = os.path.join(tmpdir, 'ace')
-            xsdir_file = os.path.join(tmpdir, 'xsdir')
-            pendf_file = os.path.join(tmpdir, 'pendf')
+            kwargs.setdefault('ace', os.path.join(tmpdir, 'ace'))
+            kwargs.setdefault('xsdir', os.path.join(tmpdir, 'xsdir'))
+            kwargs.setdefault('pendf', os.path.join(tmpdir, 'pendf'))
             kwargs['evaluation'] = evaluation
-            make_ace(filename, temperatures, ace_file, xsdir_file,
-                     pendf_file, **kwargs)
+            make_ace(filename, temperatures, **kwargs)
 
             # Create instance from ACE tables within library
-            lib = Library(ace_file)
+            lib = Library(kwargs['ace'])
             data = cls.from_ace(lib.tables[0])
             for table in lib.tables[1:]:
                 data.add_temperature_from_ace(table)
@@ -821,7 +820,7 @@ class IncidentNeutron(EqualityMixin):
 
             # Add 0K elastic scattering cross section
             if '0K' not in data.energy:
-                pendf = Evaluation(pendf_file)
+                pendf = Evaluation(kwargs['pendf'])
                 file_obj = StringIO(pendf.section[3, 2])
                 get_head_record(file_obj)
                 params, xs = get_tab1_record(file_obj)

openmc/data/thermal.py

@@ -761,15 +761,14 @@ class ThermalScattering(EqualityMixin):
         """
         with tempfile.TemporaryDirectory() as tmpdir:
             # Run NJOY to create an ACE library
-            ace_file = os.path.join(tmpdir, 'ace')
-            xsdir_file = os.path.join(tmpdir, 'xsdir')
+            kwargs.setdefault('ace', os.path.join(tmpdir, 'ace'))
+            kwargs.setdefault('xsdir', os.path.join(tmpdir, 'xsdir'))
             kwargs['evaluation'] = evaluation
             kwargs['evaluation_thermal'] = evaluation_thermal
-            make_ace_thermal(filename, filename_thermal, temperatures,
-                             ace_file, xsdir_file, **kwargs)
+            make_ace_thermal(filename, filename_thermal, temperatures, **kwargs)
 
             # Create instance from ACE tables within library
-            lib = Library(ace_file)
+            lib = Library(kwargs['ace'])
             data = cls.from_ace(lib.tables[0])
             for table in lib.tables[1:]:
                 data.add_temperature_from_ace(table)

openmc/model/funcs.py

@@ -409,10 +409,10 @@ def cylinder_from_points(p1, p2, r, **kwargs):
     kwargs['d'] = -2*dx*dy
     kwargs['e'] = -2*dy*dz
     kwargs['f'] = -2*dx*dz
-    kwargs['g'] = cy*dz - cz*dy
-    kwargs['h'] = cz*dx - cx*dz
-    kwargs['j'] = cx*dy - cy*dx
-    kwargs['k'] = -(dx*dx + dy*dy + dz*dz)*r*r
+    kwargs['g'] = 2*(cy*dz - cz*dy)
+    kwargs['h'] = 2*(cz*dx - cx*dz)
+    kwargs['j'] = 2*(cx*dy - cy*dx)
+    kwargs['k'] = cx*cx + cy*cy + cz*cz - (dx*dx + dy*dy + dz*dz)*r*r
 
     return openmc.Quadric(**kwargs)
 

tests/unit_tests/test_surface.py

@@ -362,3 +362,31 @@ def test_cylinder_from_points():
         assert p2 + 1.1*r*n in +s
         assert p1 + 0.9*r*n in -s
         assert p2 + 0.9*r*n in -s
+
+
+def test_cylinder_from_points_axis():
+    # Create axis-aligned cylinders and confirm the coefficients are as expected
+
+    # (x - 3)^2 + (y - 4)^2 = 2^2
+    # x^2 + y^2 - 6x - 8y + 21 = 0
+    s = openmc.model.cylinder_from_points((3., 4., 0.), (3., 4., 1.), 2.)
+    assert (s.a, s.b, s.c) == pytest.approx((1., 1., 0.))
+    assert (s.d, s.e, s.f) == pytest.approx((0., 0., 0.))
+    assert (s.g, s.h, s.j) == pytest.approx((-6., -8., 0.))
+    assert s.k == pytest.approx(21.)
+
+    # (y + 7)^2 + (z - 1)^2 = 3^2
+    # y^2 + z^2 + 14y - 2z + 41 = 0
+    s = openmc.model.cylinder_from_points((0., -7, 1.), (1., -7., 1.), 3.)
+    assert (s.a, s.b, s.c) == pytest.approx((0., 1., 1.))
+    assert (s.d, s.e, s.f) == pytest.approx((0., 0., 0.))
+    assert (s.g, s.h, s.j) == pytest.approx((0., 14., -2.))
+    assert s.k == 41.
+
+    # (x - 2)^2 + (z - 5)^2 = 4^2
+    # x^2 + z^2 - 4x - 10z + 13 = 0
+    s = openmc.model.cylinder_from_points((2., 0., 5.), (2., 1., 5.), 4.)
+    assert (s.a, s.b, s.c) == pytest.approx((1., 0., 1.))
+    assert (s.d, s.e, s.f) == pytest.approx((0., 0., 0.))
+    assert (s.g, s.h, s.j) == pytest.approx((-4., 0., -10.))
+    assert s.k == pytest.approx(13.)