diff --git a/test/test_new_world_grudge.py b/test/test_new_world_grudge.py
new file mode 100644
index 0000000000000000000000000000000000000000..82bfb3e2d9027e358fe2c92cb057c08026896af9
--- /dev/null
+++ b/test/test_new_world_grudge.py
@@ -0,0 +1,119 @@
+__copyright__ = """
+Copyright (C) 2021 University of Illinois Board of Trustees
+"""
+
+__license__ = """
+Permission is hereby granted, free of charge, to any person obtaining a copy
+of this software and associated documentation files (the "Software"), to deal
+in the Software without restriction, including without limitation the rights
+to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+copies of the Software, and to permit persons to whom the Software is
+furnished to do so, subject to the following conditions:
+
+The above copyright notice and this permission notice shall be included in
+all copies or substantial portions of the Software.
+
+THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
+THE SOFTWARE.
+"""
+
+import numpy as np
+import numpy.linalg as la
+
+from meshmode import _acf       # noqa: F401
+from meshmode.dof_array import flatten, thaw
+from meshmode.array_context import (  # noqa
+        pytest_generate_tests_for_pyopencl_array_context
+        as pytest_generate_tests)
+import meshmode.mesh.generation as mgen
+
+from pytools.obj_array import flat_obj_array, make_obj_array
+
+from grudge import DiscretizationCollection
+
+import pytest
+
+import logging
+
+logger = logging.getLogger(__name__)
+
+
+# {{{ mass operator trig integration
+
+@pytest.mark.parametrize("ambient_dim", [1, 2, 3])
+@pytest.mark.parametrize("quad_tag", [sym.QTAG_NONE, "OVSMP"])
+def test_mass_mat_trig(actx_factory, ambient_dim, quad_tag):
+    """Check the integral of some trig functions on an interval using the mass
+    matrix.
+    """
+    actx = actx_factory()
+
+    nelements = 17
+    order = 4
+
+    a = -4.0 * np.pi
+    b = +9.0 * np.pi
+    true_integral = 13*np.pi/2 * (b - a)**(ambient_dim - 1)
+
+    from meshmode.discretization.poly_element import \
+        QuadratureSimplexGroupFactory
+
+    dd_quad = sym.DOFDesc(sym.DTAG_VOLUME_ALL, quad_tag)
+
+    if quad_tag is sym.QTAG_NONE:
+        quad_tag_to_group_factory = {}
+    else:
+        quad_tag_to_group_factory = {
+            quad_tag: QuadratureSimplexGroupFactory(order=2*order)
+        }
+
+    mesh = mgen.generate_regular_rect_mesh(
+        a=(a,)*ambient_dim, b=(b,)*ambient_dim,
+        n=(nelements,)*ambient_dim,
+        order=1
+    )
+    dcoll = DiscretizationCollection(
+        actx, mesh, order=order,
+        quad_tag_to_group_factory=quad_tag_to_group_factory
+    )
+
+    def f(x):
+        return actx.np.sin(x[0])**2
+
+    def ones(x):
+        return actx.np.ones(x.shape)
+
+    volm_disc = dcoll.disc_from_dd(sym.DD_VOLUME)
+    x_volm = thaw(actx, volm_disc.nodes())
+    f_volm = f(x_vol)
+    ones_volm = ones(x_vol)
+
+    quad_disc = dcoll.disc_from_dd(dd_quad)
+    x_quad = thaw(actx, quad_disc.nodes())
+    f_quad = f(x_quad)
+    ones_quad = ones(x_quad)
+
+    mass_op = bind(discr, sym.MassOperator(dd_quad, sym.DD_VOLUME)(sym_f))
+
+    num_integral_1 = np.dot(ones_volm, actx.to_numpy(flatten(mass_op(f=f_quad))))
+    err_1 = abs(num_integral_1 - true_integral)
+    assert err_1 < 1e-9, err_1
+
+    num_integral_2 = np.dot(f_volm, actx.to_numpy(flatten(mass_op(f=ones_quad))))
+    err_2 = abs(num_integral_2 - true_integral)
+    assert err_2 < 1.0e-9, err_2
+
+    if quad_tag is sym.QTAG_NONE:
+        # NOTE: `integral` always makes a square mass matrix and
+        # `QuadratureSimplexGroupFactory` does not have a `mass_matrix` method.
+        num_integral_3 = bind(discr,
+                sym.integral(sym_f, dd=dd_quad))(f=f_quad)
+        err_3 = abs(num_integral_3 - true_integral)
+        assert err_3 < 5.0e-10, err_3
+
+# }}}
\ No newline at end of file