diff --git a/tests/regression_tests/torus/__init__.py b/tests/regression_tests/torus/__init__.py
new file mode 100644
index 0000000000..e69de29bb2
diff --git a/tests/regression_tests/torus/inputs_true.dat b/tests/regression_tests/torus/inputs_true.dat
new file mode 100644
index 0000000000..0045d7e466
--- /dev/null
+++ b/tests/regression_tests/torus/inputs_true.dat
@@ -0,0 +1,34 @@
+<?xml version='1.0' encoding='utf-8'?>
+<geometry>
+  <cell id="1" material="1" region="-2" universe="1" />
+  <cell id="2" material="1" region="-3" universe="1" />
+  <cell id="3" material="1" region="-1" universe="1" />
+  <cell id="4" material="2" region="4 -5 6 -7 8 -9 2 3 1" universe="1" />
+  <surface coeffs="0.0 0.0 0.0 3 1.5 1" id="1" type="z-torus" />
+  <surface coeffs="6 0.0 0.0 3 1.5 1" id="2" type="x-torus" />
+  <surface coeffs="6 0.0 0.0 6 1 0.75" id="3" type="y-torus" />
+  <surface boundary="vacuum" coeffs="-5" id="4" type="x-plane" />
+  <surface boundary="vacuum" coeffs="14" id="5" type="x-plane" />
+  <surface boundary="vacuum" coeffs="-5" id="6" type="y-plane" />
+  <surface boundary="vacuum" coeffs="5" id="7" type="y-plane" />
+  <surface boundary="vacuum" coeffs="-8" id="8" type="z-plane" />
+  <surface boundary="vacuum" coeffs="8" id="9" type="z-plane" />
+</geometry>
+<?xml version='1.0' encoding='utf-8'?>
+<materials>
+  <material depletable="true" id="1">
+    <density units="g/cm3" value="12.0" />
+    <nuclide ao="1.0" name="U235" />
+  </material>
+  <material id="2">
+    <density units="g/cm3" value="1.0" />
+    <nuclide ao="1.0" name="H1" />
+  </material>
+</materials>
+<?xml version='1.0' encoding='utf-8'?>
+<settings>
+  <run_mode>eigenvalue</run_mode>
+  <particles>1000</particles>
+  <batches>10</batches>
+  <inactive>5</inactive>
+</settings>
diff --git a/tests/regression_tests/torus/results_true.dat b/tests/regression_tests/torus/results_true.dat
new file mode 100644
index 0000000000..42444ad09a
--- /dev/null
+++ b/tests/regression_tests/torus/results_true.dat
@@ -0,0 +1,2 @@
+k-combined:
+7.628424E-01 2.557300E-02
diff --git a/tests/regression_tests/torus/test.py b/tests/regression_tests/torus/test.py
new file mode 100644
index 0000000000..8970a5cc77
--- /dev/null
+++ b/tests/regression_tests/torus/test.py
@@ -0,0 +1,39 @@
+import openmc
+import pytest
+
+from tests.testing_harness import PyAPITestHarness
+
+
+@pytest.fixture
+def model():
+    model = openmc.Model()
+    fuel = openmc.Material()
+    fuel.set_density('g/cm3', 12.0)
+    fuel.add_nuclide('U235', 1.0)
+    al = openmc.Material()
+    al.set_density('g/cm3', 1.0)
+    al.add_nuclide('H1', 1.0)
+    model.materials.extend([fuel, al])
+
+    # 🍩🍩🍩
+    zt = openmc.ZTorus(a=3, b=1.5, c=1)
+    xt = openmc.XTorus(x0=6, a=3, b=1.5, c=1)
+    yt = openmc.YTorus(x0=6, a=6, b=1, c=0.75)
+    box = openmc.model.RectangularParallelepiped(-5, 14, -5, 5, -8, 8,
+                                                 boundary_type='vacuum')
+
+    xt_cell = openmc.Cell(fill=fuel, region=-xt)
+    yt_cell = openmc.Cell(fill=fuel, region=-yt)
+    zt_cell = openmc.Cell(fill=fuel, region=-zt)
+    outer_cell = openmc.Cell(fill=al, region=-box & +xt & +yt & +zt)
+    model.geometry = openmc.Geometry([xt_cell, yt_cell, zt_cell, outer_cell])
+
+    model.settings.particles = 1000
+    model.settings.batches = 10
+    model.settings.inactive = 5
+    return model
+
+
+def test_torus(model):
+    harness = PyAPITestHarness('statepoint.10.h5', model)
+    harness.main()
diff --git a/tests/unit_tests/test_surface.py b/tests/unit_tests/test_surface.py
index e7d37c4e91..c24d6e3521 100644
--- a/tests/unit_tests/test_surface.py
+++ b/tests/unit_tests/test_surface.py
@@ -605,3 +605,127 @@ def test_cylinder_from_points_axis():
     assert (d, e, f) == pytest.approx((0., 0., 0.))
     assert (g, h, j) == pytest.approx((-4., 0., -10.))
     assert k == pytest.approx(13.)
+
+
+def torus_common(center, R, r1, r2, cls):
+    x, y, z = center
+    s = cls(x0=x, y0=y, z0=z, a=R, b=r1, c=r2)
+    assert s.x0 == x
+    assert s.y0 == y
+    assert s.z0 == z
+    assert s.a == R
+    assert s.b == r1
+    assert s.c == r2
+
+    # evaluate method
+    assert s.evaluate((x, y, z)) > 0.0
+
+    # translate method
+    trans = np.array([1.0, 1.5, -2.0])
+    st = s.translate(trans)
+    assert st.x0 == s.x0 + trans[0]
+    assert st.y0 == s.y0 + trans[1]
+    assert st.z0 == s.z0 + trans[2]
+    assert st.a == s.a
+    assert st.b == s.b
+    assert st.c == s.c
+
+    # can't rotate at present
+    with pytest.raises(NotImplementedError):
+        s.rotate((0., 1., 0.))
+
+    # Check bounding box
+    ll, ur = (+s).bounding_box
+    assert np.all(np.isinf(ll))
+    assert np.all(np.isinf(ur))
+
+    ll, ur = (-s).bounding_box
+    llt, urt = (-st).bounding_box
+    np.testing.assert_allclose(ll + trans, llt)
+    np.testing.assert_allclose(ur + trans, urt)
+
+    # Make sure repr works
+    repr(s)
+
+    return s
+
+
+def test_xtorus():
+    x, y, z = 2, -4, 5
+    R, r1, r2 = 3, 1.5, 1
+    s = torus_common((x, y, z), R, r1, r2, openmc.XTorus)
+
+    # evaluate method (points inside torus)
+    assert s.evaluate((x, y + R, z)) < 0.0
+    assert s.evaluate((x, y - R, z)) < 0.0
+    assert s.evaluate((x, y, z + R)) < 0.0
+    assert s.evaluate((x, y, z - R)) < 0.0
+
+    # evaluate method (points on torus)
+    assert s.evaluate((x, y + R + r2, z)) == pytest.approx(0.0)
+    assert s.evaluate((x, y - R - r2, z)) == pytest.approx(0.0)
+    assert s.evaluate((x, y, z + R - r2)) == pytest.approx(0.0)
+    assert s.evaluate((x, y, z - R + r2)) == pytest.approx(0.0)
+    assert s.evaluate((x + r1, y + R, z)) == pytest.approx(0.0)
+
+    # evaluate method (points outside torus)
+    assert s.evaluate((x, y + R, z + R)) > 0.0
+    assert s.evaluate((x, y - R, z + R)) > 0.0
+    assert s.evaluate((x, y + R + r2 + 0.01, z)) > 0.0
+    assert s.evaluate((x, y, z + R + r2 + 0.01)) > 0.0
+    assert s.evaluate((x + r1 + 0.01, y, z + R)) > 0.0
+    assert s.evaluate((x + r1 + 0.01, y + R, z)) > 0.0
+
+
+def test_ytorus():
+    x, y, z = 2, -4, 5
+    R, r1, r2 = 3, 1.5, 1
+    s = torus_common((x, y, z), R, r1, r2, openmc.YTorus)
+
+    # evaluate method (points inside torus)
+    assert s.evaluate((x + R, y, z)) < 0.0
+    assert s.evaluate((x - R, y, z)) < 0.0
+    assert s.evaluate((x, y, z + R)) < 0.0
+    assert s.evaluate((x, y, z - R)) < 0.0
+
+    # evaluate method (points on torus)
+    assert s.evaluate((x + R + r2, y, z)) == pytest.approx(0.0)
+    assert s.evaluate((x - R - r2, y, z)) == pytest.approx(0.0)
+    assert s.evaluate((x, y, z + R - r2)) == pytest.approx(0.0)
+    assert s.evaluate((x, y, z - R + r2)) == pytest.approx(0.0)
+    assert s.evaluate((x + R, y + r1, z)) == pytest.approx(0.0)
+
+    # evaluate method (points outside torus)
+    assert s.evaluate((x + R, y, z + R)) > 0.0
+    assert s.evaluate((x - R, y, z + R)) > 0.0
+    assert s.evaluate((x + R + r2 + 0.01, y, z)) > 0.0
+    assert s.evaluate((x, y, z + R + r2 + 0.01)) > 0.0
+    assert s.evaluate((x, y + r1 + 0.01, z + R)) > 0.0
+    assert s.evaluate((x + R, y + r1 + 0.01, z)) > 0.0
+
+
+def test_ztorus():
+    x, y, z = 2, -4, 5
+    R, r1, r2 = 3, 1.5, 1
+    s = torus_common((x, y, z), R, r1, r2, openmc.ZTorus)
+
+    # evaluate method (points inside torus)
+    assert s.evaluate((x, y + R, z)) < 0.0
+    assert s.evaluate((x, y - R, z)) < 0.0
+    assert s.evaluate((x + R, y, z)) < 0.0
+    assert s.evaluate((x - R, y, z)) < 0.0
+
+    # evaluate method (points on torus)
+    assert s.evaluate((x, y + R + r2, z)) == pytest.approx(0.0)
+    assert s.evaluate((x, y - R - r2, z)) == pytest.approx(0.0)
+    assert s.evaluate((x + R - r2, y, z)) == pytest.approx(0.0)
+    assert s.evaluate((x - R + r2, y, z)) == pytest.approx(0.0)
+    assert s.evaluate((x, y + R, z + r1)) == pytest.approx(0.0)
+
+    # evaluate method (points outside torus)
+    assert s.evaluate((x + R, y + R, z)) > 0.0
+    assert s.evaluate((x + R, y - R, z)) > 0.0
+    assert s.evaluate((x, y + R + r2 + 0.01, z)) > 0.0
+    assert s.evaluate((x + R + r2 + 0.01, y, z)) > 0.0
+    assert s.evaluate((x + R, y, z + r1 + 0.01)) > 0.0
+    assert s.evaluate((x, y + R, z + r1 + 0.01)) > 0.0
diff --git a/tests/unit_tests/test_torus.py b/tests/unit_tests/test_torus.py
new file mode 100644
index 0000000000..925949afa7
--- /dev/null
+++ b/tests/unit_tests/test_torus.py
@@ -0,0 +1,45 @@
+from itertools import combinations
+from random import uniform
+import openmc
+
+
+def get_torus_keff(cls, center=(0, 0, 0)):
+    model = openmc.Model()
+    mat = openmc.Material()
+    mat.add_nuclide('U235', 1.0)
+    mat.set_density('g/cm3', 10.0)
+
+    x, y, z = center
+    torus = cls(x0=x, y0=y, z0=z, a=3, b=2, c=2)
+    sphere = openmc.Sphere(x0=x, y0=y, z0=z, r=5.5, boundary_type="vacuum")
+    torus_cell = openmc.Cell(fill=mat, region=-torus)
+    outer_cell = openmc.Cell(region=+torus & -sphere)
+    model.geometry = openmc.Geometry([torus_cell, outer_cell])
+
+    model.settings.source = openmc.Source(space=openmc.stats.Point(center))
+    model.settings.batches = 10
+    model.settings.inactive = 5
+    model.settings.particles = 1000
+
+    sp_path = model.run()
+    with openmc.StatePoint(sp_path) as sp:
+        return sp.k_combined
+
+
+def test_torus_keff(run_in_tmpdir):
+    random_point = lambda: (uniform(-5, 5), uniform(-5, 5), uniform(-5, 5))
+    keffs = [
+        get_torus_keff(openmc.XTorus),
+        get_torus_keff(openmc.XTorus, random_point()),
+        get_torus_keff(openmc.YTorus),
+        get_torus_keff(openmc.YTorus, random_point()),
+        get_torus_keff(openmc.ZTorus),
+        get_torus_keff(openmc.ZTorus, random_point())
+    ]
+
+    # For each combination of keff values, their difference should be within
+    # uncertainty (3 std dev)
+    for k1, k2 in combinations(keffs, 2):
+        print(k1, k2)
+        diff = k1 - k2
+        assert abs(diff.n) < 3*diff.s