diff --git a/test/mesh_data.py b/test/mesh_data.py
index 0671cd87d0beac099a3c50625d0cc68d878c60e7..0ccc369a1347c01517b3a8a30089558aa4f4aaf7 100644
--- a/test/mesh_data.py
+++ b/test/mesh_data.py
@@ -61,8 +61,8 @@ class SphereMeshBuilder(MeshBuilder):
     radius = 1.0
 
     def get_mesh(self, resolution, mesh_order):
-        from meshmode.mesh.generation import generate_icosphere
-        return generate_icosphere(self.radius, order=mesh_order,
+        from meshmode.mesh.generation import generate_sphere
+        return generate_sphere(self.radius, order=mesh_order,
                 uniform_refinement_rounds=resolution)
 
 
@@ -70,41 +70,18 @@ class SpheroidMeshBuilder(MeshBuilder):
     ambient_dim = 3
 
     mesh_order = 4
-    resolutions = [1.0, 0.11, 0.05]
-    # resolutions = [1.0, 0.11, 0.05, 0.03, 0.015]
-
-    @property
-    def radius(self):
-        return 0.25
-
-    @property
-    def diameter(self):
-        return 2.0 * self.radius
+    resolutions = [0, 1, 2, 3]
 
-    @property
-    def aspect_ratio(self):
-        return 2.0
+    radius = 1.0
+    aspect_ratio = 2.0
 
     def get_mesh(self, resolution, mesh_order):
-        from meshmode.mesh.io import ScriptSource
-        source = ScriptSource("""
-            SetFactory("OpenCASCADE");
-            Sphere(1) = {{0, 0, 0, {r}}};
-            Dilate {{ {{0, 0, 0}}, {{ {r}, {r}, {rr} }} }} {{ Volume{{1}}; }}
-            """.format(r=self.diameter, rr=self.aspect_ratio * self.diameter),
-            "geo"
-        )
-
-        from meshmode.mesh.io import generate_gmsh
-        mesh = generate_gmsh(source, 2, order=mesh_order,
-                other_options=[
-                    "-optimize_ho",
-                    "-string", "Mesh.CharacteristicLengthMax = %g;" % resolution
-                    ],
-                target_unit="MM")
-
-        from meshmode.mesh.processing import perform_flips
-        return perform_flips(mesh, np.ones(mesh.nelements))
+        from meshmode.mesh.generation import generate_sphere
+        mesh = generate_sphere(self.radius, order=mesh_order,
+                uniform_refinement_rounds=resolution)
+
+        from meshmode.mesh.processing import affine_map
+        return affine_map(mesh, A=np.diag([1.0, 1.0, self.aspect_ratio]))
 
 
 class BoxMeshBuilder(MeshBuilder):
diff --git a/test/test_grudge.py b/test/test_grudge.py
index 65f31336943917a865752d43c96e6cb463015a5e..9de657e1c619acdbd06c47fe6b7e39573c7c0c79 100644
--- a/test/test_grudge.py
+++ b/test/test_grudge.py
@@ -33,6 +33,7 @@ pytest_generate_tests = pytest_generate_tests_for_array_contexts(
 
 from arraycontext.container.traversal import thaw
 
+from meshmode import _acf           # noqa: F401
 from meshmode.dof_array import flat_norm
 import meshmode.mesh.generation as mgen
 
@@ -509,7 +510,7 @@ def test_2d_gauss_theorem(actx_factory):
     assert abs(int_1 - int_2) < 1e-13
 
 
-@pytest.mark.parametrize("mesh_name", ["2-1-ellipse", "spheroid"])
+@pytest.mark.parametrize("mesh_name", ["2-1-ellipse", "2-1-spheroid"])
 def test_surface_divergence_theorem(actx_factory, mesh_name, visualize=False):
     r"""Check the surface divergence theorem.
 
@@ -531,9 +532,9 @@ def test_surface_divergence_theorem(actx_factory, mesh_name, visualize=False):
     if mesh_name == "2-1-ellipse":
         from mesh_data import EllipseMeshBuilder
         builder = EllipseMeshBuilder(radius=3.1, aspect_ratio=2.0)
-    elif mesh_name == "spheroid":
+    elif mesh_name == "2-1-spheroid":
         from mesh_data import SpheroidMeshBuilder
-        builder = SpheroidMeshBuilder()
+        builder = SpheroidMeshBuilder(radius=3.1, aspect_ratio=2.0)
     elif mesh_name == "circle":
         from mesh_data import EllipseMeshBuilder
         builder = EllipseMeshBuilder(radius=1.0, aspect_ratio=1.0)
@@ -627,17 +628,17 @@ def test_surface_divergence_theorem(actx_factory, mesh_name, visualize=False):
         op_global = op.nodal_sum(dcoll, dd, stiff - (stiff_t + kterm))
         op_local = op.elementwise_sum(dcoll, dd, stiff - (stiff_t + kterm + flux))
 
-        err_global = abs(op_global)
-        err_local = op.norm(dcoll, op_local, np.inf)
-        logger.info("errors: global %.5e local %.5e", err_global, err_local)
-
         # compute max element size
         from grudge.dt_utils import h_max_from_volume
 
-        h_max = h_max_from_volume(dcoll)
+        h_max = actx.to_numpy(h_max_from_volume(dcoll))
+        err_global = actx.to_numpy(abs(op_global))
+        err_local = actx.to_numpy(op.norm(dcoll, op_local, np.inf))
+        logger.info("errors: h_max %.5e global %.5e local %.5e",
+                h_max, err_global, err_local)
 
-        eoc_global.add_data_point(actx.to_numpy(h_max), actx.to_numpy(err_global))
-        eoc_local.add_data_point(actx.to_numpy(h_max), actx.to_numpy(err_local))
+        eoc_global.add_data_point(h_max, err_global)
+        eoc_local.add_data_point(h_max, err_local)
 
         if visualize:
             from grudge.shortcuts import make_visualizer
@@ -651,8 +652,8 @@ def test_surface_divergence_theorem(actx_factory, mesh_name, visualize=False):
     # }}}
 
     order = min(builder.order, builder.mesh_order) - 0.5
-    logger.info("\n%s", str(eoc_global))
-    logger.info("\n%s", str(eoc_local))
+    logger.info("eoc_global:\n%s", str(eoc_global))
+    logger.info("eoc_local:\n%s", str(eoc_local))
 
     assert eoc_global.max_error() < 1.0e-12 \
             or eoc_global.order_estimate() > order - 0.5