diff --git a/examples/lump-mpi.py b/examples/lump-mpi.py
index 3df2bef1e..5f570f1ea 100644
--- a/examples/lump-mpi.py
+++ b/examples/lump-mpi.py
@@ -29,6 +29,7 @@
 import pyopencl.tools as cl_tools
 from functools import partial
 
+from arraycontext import flatten
 from meshmode.array_context import (
     PyOpenCLArrayContext,
     PytatoPyOpenCLArrayContext
@@ -43,7 +44,8 @@
 from mirgecom.euler import euler_operator
 from mirgecom.simutil import (
     get_sim_timestep,
-    generate_and_distribute_mesh
+    generate_and_distribute_mesh,
+    componentwise_norms
 )
 from mirgecom.io import make_init_message
 from mirgecom.mpi import mpi_entry_point
@@ -247,10 +249,10 @@ def my_health_check(dv, state, exact):
             health_error = True
             logger.info(f"{rank=}: Invalid pressure data found.")
 
-        from mirgecom.simutil import compare_fluid_solutions
-        component_errors = compare_fluid_solutions(discr, state, exact)
+        component_errors = componentwise_norms(discr, state - exact)
+        max_error = actx.to_numpy(actx.np.max(component_errors))
         exittol = .09
-        if max(component_errors) > exittol:
+        if max_error > exittol:
             health_error = True
             if rank == 0:
                 logger.info("Solution diverged from exact soln.")
@@ -296,11 +298,11 @@ def my_pre_step(step, t, dt, state):
             if do_status:
                 if exact is None:
                     exact = initializer(x_vec=nodes, eos=eos, time=t)
-                from mirgecom.simutil import compare_fluid_solutions
-                component_errors = compare_fluid_solutions(discr, state, exact)
+                flat_component_errors = actx.to_numpy(flatten(
+                    componentwise_norms(discr, state - exact), actx))
                 status_msg = (
                     "------- errors="
-                    + ", ".join("%.3g" % en for en in component_errors))
+                    + ", ".join("%.3g" % en for en in flat_component_errors))
                 if rank == 0:
                     logger.info(status_msg)
 
diff --git a/examples/mixture-mpi.py b/examples/mixture-mpi.py
index a0905c017..0a5540575 100644
--- a/examples/mixture-mpi.py
+++ b/examples/mixture-mpi.py
@@ -29,6 +29,7 @@
 import pyopencl.tools as cl_tools
 from functools import partial
 
+from arraycontext import flatten
 from meshmode.array_context import (
     PyOpenCLArrayContext,
     PytatoPyOpenCLArrayContext
@@ -43,7 +44,8 @@
 from mirgecom.euler import euler_operator
 from mirgecom.simutil import (
     get_sim_timestep,
-    generate_and_distribute_mesh
+    generate_and_distribute_mesh,
+    componentwise_norms
 )
 from mirgecom.io import make_init_message
 from mirgecom.mpi import mpi_entry_point
@@ -225,9 +227,10 @@ def main(ctx_factory=cl.create_some_context, use_logmgr=True,
         logger.info(init_message)
 
     def my_write_status(component_errors):
+        flat_component_errors = actx.to_numpy(flatten(component_errors, actx))
         status_msg = (
             "------- errors="
-            + ", ".join("%.3g" % en for en in component_errors))
+            + ", ".join("%.3g" % en for en in flat_component_errors))
         if rank == 0:
             logger.info(status_msg)
 
@@ -270,8 +273,9 @@ def my_health_check(dv, component_errors):
             health_error = True
             logger.info(f"{rank=}: Invalid pressure data found.")
 
+        max_error = actx.to_numpy(actx.np.max(component_errors))
         exittol = .09
-        if max(component_errors) > exittol:
+        if max_error > exittol:
             health_error = True
             if rank == 0:
                 logger.info("Solution diverged from exact soln.")
@@ -296,8 +300,7 @@ def my_pre_step(step, t, dt, state):
             if do_health:
                 dv = eos.dependent_vars(state)
                 exact = initializer(x_vec=nodes, eos=eos, time=t)
-                from mirgecom.simutil import compare_fluid_solutions
-                component_errors = compare_fluid_solutions(discr, state, exact)
+                component_errors = componentwise_norms(discr, state - exact)
                 health_errors = global_reduce(
                     my_health_check(dv, component_errors), op="lor")
                 if health_errors:
@@ -321,8 +324,7 @@ def my_pre_step(step, t, dt, state):
                 if component_errors is None:
                     if exact is None:
                         exact = initializer(x_vec=nodes, eos=eos, time=t)
-                    from mirgecom.simutil import compare_fluid_solutions
-                    component_errors = compare_fluid_solutions(discr, state, exact)
+                    component_errors = componentwise_norms(discr, state - exact)
                 my_write_status(component_errors)
 
         except MyRuntimeError:
diff --git a/examples/scalar-lump-mpi.py b/examples/scalar-lump-mpi.py
index d0861e3e3..51423d442 100644
--- a/examples/scalar-lump-mpi.py
+++ b/examples/scalar-lump-mpi.py
@@ -30,6 +30,7 @@
 from functools import partial
 from pytools.obj_array import make_obj_array
 
+from arraycontext import flatten
 from meshmode.array_context import (
     PyOpenCLArrayContext,
     PytatoPyOpenCLArrayContext
@@ -44,7 +45,8 @@
 from mirgecom.euler import euler_operator
 from mirgecom.simutil import (
     get_sim_timestep,
-    generate_and_distribute_mesh
+    generate_and_distribute_mesh,
+    componentwise_norms
 )
 from mirgecom.io import make_init_message
 from mirgecom.mpi import mpi_entry_point
@@ -215,10 +217,11 @@ def main(ctx_factory=cl.create_some_context, use_logmgr=True,
         logger.info(init_message)
 
     def my_write_status(component_errors):
+        flat_component_errors = actx.to_numpy(flatten(component_errors, actx))
         if rank == 0:
             logger.info(
                 "------- errors="
-                + ", ".join("%.3g" % en for en in component_errors))
+                + ", ".join("%.3g" % en for en in flat_component_errors))
 
     def my_write_viz(step, t, state, dv=None, exact=None, resid=None):
         if dv is None:
@@ -258,8 +261,9 @@ def my_health_check(pressure, component_errors):
             health_error = True
             logger.info(f"{rank=}: Invalid pressure data found.")
 
+        max_error = actx.to_numpy(actx.np.max(component_errors))
         exittol = .09
-        if max(component_errors) > exittol:
+        if max_error > exittol:
             health_error = True
             if rank == 0:
                 logger.info("Solution diverged from exact soln.")
@@ -284,8 +288,7 @@ def my_pre_step(step, t, dt, state):
             if do_health:
                 dv = eos.dependent_vars(state)
                 exact = initializer(x_vec=nodes, eos=eos, time=t)
-                from mirgecom.simutil import compare_fluid_solutions
-                component_errors = compare_fluid_solutions(discr, state, exact)
+                component_errors = componentwise_norms(discr, state - exact)
                 health_errors = global_reduce(
                     my_health_check(dv.pressure, component_errors), op="lor")
                 if health_errors:
@@ -309,8 +312,7 @@ def my_pre_step(step, t, dt, state):
                 if component_errors is None:
                     if exact is None:
                         exact = initializer(x_vec=nodes, eos=eos, time=t)
-                    from mirgecom.simutil import compare_fluid_solutions
-                    component_errors = compare_fluid_solutions(discr, state, exact)
+                    component_errors = componentwise_norms(discr, state - exact)
                 my_write_status(component_errors)
 
         except MyRuntimeError:
diff --git a/examples/sod-mpi.py b/examples/sod-mpi.py
index 755ac2212..8ff8828c9 100644
--- a/examples/sod-mpi.py
+++ b/examples/sod-mpi.py
@@ -29,6 +29,7 @@
 import pyopencl.tools as cl_tools
 from functools import partial
 
+from arraycontext import flatten
 from meshmode.array_context import (
     PyOpenCLArrayContext,
     PytatoPyOpenCLArrayContext
@@ -43,7 +44,8 @@
 from mirgecom.euler import euler_operator
 from mirgecom.simutil import (
     get_sim_timestep,
-    generate_and_distribute_mesh
+    generate_and_distribute_mesh,
+    componentwise_norms
 )
 from mirgecom.io import make_init_message
 from mirgecom.mpi import mpi_entry_point
@@ -204,10 +206,11 @@ def main(ctx_factory=cl.create_some_context, use_logmgr=True,
         logger.info(init_message)
 
     def my_write_status(component_errors):
+        flat_component_errors = actx.to_numpy(flatten(component_errors, actx))
         if rank == 0:
             logger.info(
                 "------- errors="
-                + ", ".join("%.3g" % en for en in component_errors)
+                + ", ".join("%.3g" % en for en in flat_component_errors)
             )
 
     def my_write_viz(step, t, state, dv=None, exact=None, resid=None):
@@ -248,8 +251,9 @@ def my_health_check(pressure, component_errors):
             health_error = True
             logger.info(f"{rank=}: Invalid pressure data found.")
 
+        max_error = actx.to_numpy(actx.np.max(component_errors))
         exittol = .09
-        if max(component_errors) > exittol:
+        if max_error > exittol:
             health_error = True
             if rank == 0:
                 logger.info("Solution diverged from exact soln.")
@@ -274,8 +278,7 @@ def my_pre_step(step, t, dt, state):
             if do_health:
                 dv = eos.dependent_vars(state)
                 exact = initializer(x_vec=nodes, eos=eos, time=t)
-                from mirgecom.simutil import compare_fluid_solutions
-                component_errors = compare_fluid_solutions(discr, state, exact)
+                component_errors = componentwise_norms(discr, state - exact)
                 health_errors = global_reduce(
                     my_health_check(dv.pressure, component_errors), op="lor")
                 if health_errors:
@@ -299,9 +302,7 @@ def my_pre_step(step, t, dt, state):
                 if component_errors is None:
                     if exact is None:
                         exact = initializer(x_vec=nodes, eos=eos, time=t)
-                        from mirgecom.simutil import compare_fluid_solutions
-                        component_errors = \
-                            compare_fluid_solutions(discr, state, exact)
+                    component_errors = componentwise_norms(discr, state - exact)
                 my_write_status(component_errors)
 
         except MyRuntimeError:
diff --git a/examples/vortex-mpi.py b/examples/vortex-mpi.py
index fc5071029..e91f73109 100644
--- a/examples/vortex-mpi.py
+++ b/examples/vortex-mpi.py
@@ -29,6 +29,7 @@
 import pyopencl.tools as cl_tools
 from functools import partial
 
+from arraycontext import flatten
 from meshmode.array_context import (
     PyOpenCLArrayContext,
     PytatoPyOpenCLArrayContext
@@ -43,7 +44,8 @@
 from mirgecom.simutil import (
     get_sim_timestep,
     generate_and_distribute_mesh,
-    check_step
+    check_step,
+    componentwise_norms
 )
 from mirgecom.io import make_init_message
 from mirgecom.mpi import mpi_entry_point
@@ -232,10 +234,11 @@ def my_write_status(state, component_errors, cfl=None):
                         discr, "vol",
                         get_inviscid_cfl(discr, eos, current_dt, cv=state)))[()]
         if rank == 0:
+            flat_component_errors = actx.to_numpy(flatten(component_errors, actx))
             logger.info(
                 f"------ {cfl=}\n"
                 "------- errors="
-                + ", ".join("%.3g" % en for en in component_errors))
+                + ", ".join("%.3g" % en for en in flat_component_errors))
 
     def my_write_viz(step, t, state, dv=None, exact=None, resid=None):
         if dv is None:
@@ -275,8 +278,9 @@ def my_health_check(pressure, component_errors):
             health_error = True
             logger.info(f"{rank=}: Invalid pressure data found.")
 
+        max_error = actx.to_numpy(actx.np.max(component_errors))
         exittol = .1
-        if max(component_errors) > exittol:
+        if max_error > exittol:
             health_error = True
             if rank == 0:
                 logger.info("Solution diverged from exact soln.")
@@ -300,8 +304,7 @@ def my_pre_step(step, t, dt, state):
             if do_health:
                 dv = eos.dependent_vars(state)
                 exact = initializer(x_vec=nodes, eos=eos, time=t)
-                from mirgecom.simutil import compare_fluid_solutions
-                component_errors = compare_fluid_solutions(discr, state, exact)
+                component_errors = componentwise_norms(discr, state - exact)
                 health_errors = global_reduce(
                     my_health_check(dv.pressure, component_errors), op="lor")
                 if health_errors:
@@ -316,8 +319,7 @@ def my_pre_step(step, t, dt, state):
                 if component_errors is None:
                     if exact is None:
                         exact = initializer(x_vec=nodes, eos=eos, time=t)
-                    from mirgecom.simutil import compare_fluid_solutions
-                    component_errors = compare_fluid_solutions(discr, state, exact)
+                    component_errors = componentwise_norms(discr, state - exact)
                 my_write_status(state, component_errors)
 
             if do_viz:
diff --git a/mirgecom/filter.py b/mirgecom/filter.py
index da09d9526..229c5bdf2 100644
--- a/mirgecom/filter.py
+++ b/mirgecom/filter.py
@@ -48,10 +48,7 @@
 import numpy as np
 import grudge.dof_desc as dof_desc
 
-from arraycontext import map_array_container
-
-from functools import partial
-
+from arraycontext import multimapped_over_array_containers
 from meshmode.dof_array import DOFArray
 
 from pytools import keyed_memoize_in
@@ -167,6 +164,7 @@ def apply_spectral_filter(actx, modal_field, discr, cutoff,
     )
 
 
+@multimapped_over_array_containers(leaf_class=DOFArray)
 def filter_modally(dcoll, dd, cutoff, mode_resp_func, field):
     """Stand-alone procedural interface to spectral filtering.
 
@@ -200,11 +198,6 @@ def filter_modally(dcoll, dd, cutoff, mode_resp_func, field):
     result: :class:`mirgecom.fluid.ConservedVars`
         An array container containing the filtered field(s).
     """
-    if not isinstance(field, DOFArray):
-        return map_array_container(
-            partial(filter_modally, dcoll, dd, cutoff, mode_resp_func), field
-        )
-
     actx = field.array_context
     dd = dof_desc.as_dofdesc(dd)
     dd_modal = dof_desc.DD_VOLUME_MODAL
diff --git a/mirgecom/simutil.py b/mirgecom/simutil.py
index f310b7a7f..0e07058dc 100644
--- a/mirgecom/simutil.py
+++ b/mirgecom/simutil.py
@@ -13,7 +13,6 @@
 
 .. autofunction:: compare_fluid_solutions
 .. autofunction:: componentwise_norms
-.. autofunction:: max_component_norm
 .. autofunction:: check_naninf_local
 .. autofunction:: check_range_local
 
@@ -50,9 +49,7 @@
 import numpy as np
 import grudge.op as op
 
-from arraycontext import map_array_container, flatten
-
-from functools import partial
+from arraycontext import multimapped_over_array_containers
 
 from meshmode.dof_array import DOFArray
 
@@ -294,42 +291,33 @@ def check_naninf_local(discr, dd, field):
 def compare_fluid_solutions(discr, red_state, blue_state):
     """Return inf norm of (*red_state* - *blue_state*) for each component.
 
+    Deprecated. Do not use in new code.
+
     .. note::
         This is a collective routine and must be called by all MPI ranks.
     """
+    from warnings import warn
+    warn("compare_fluid_solutions is deprecated and will disappear in Q3 2022. "
+         "Use componentwise_norms instead.", DeprecationWarning, stacklevel=2)
+
     actx = red_state.array_context
     resid = red_state - blue_state
+    from arraycontext import flatten
     resid_errs = actx.to_numpy(
         flatten(componentwise_norms(discr, resid, order=np.inf), actx))
 
     return resid_errs.tolist()
 
 
+# FIXME: Add componentwise norm functionality to grudge?
+@multimapped_over_array_containers(leaf_class=DOFArray)
 def componentwise_norms(discr, fields, order=np.inf):
     """Return the *order*-norm for each component of *fields*.
 
     .. note::
         This is a collective routine and must be called by all MPI ranks.
     """
-    if not isinstance(fields, DOFArray):
-        return map_array_container(
-            partial(componentwise_norms, discr, order=order), fields)
-    if len(fields) > 0:
-        return discr.norm(fields, order)
-    else:
-        # FIXME: This work-around for #575 can go away after #569
-        return 0
-
-
-def max_component_norm(discr, fields, order=np.inf):
-    """Return the max *order*-norm over the components of *fields*.
-
-    .. note::
-        This is a collective routine and must be called by all MPI ranks.
-    """
-    actx = fields.array_context
-    return max(actx.to_numpy(flatten(
-        componentwise_norms(discr, fields, order), actx)))
+    return discr.norm(fields, order)
 
 
 def generate_and_distribute_mesh(comm, generate_mesh):
diff --git a/requirements.txt b/requirements.txt
index 4dc2b033d..006a4015c 100644
--- a/requirements.txt
+++ b/requirements.txt
@@ -16,9 +16,9 @@ pyyaml
 --editable git+https://github.com/inducer/dagrt.git#egg=dagrt
 --editable git+https://github.com/inducer/leap.git#egg=leap
 --editable git+https://github.com/inducer/modepy.git#egg=modepy
---editable git+https://github.com/inducer/arraycontext.git#egg=arraycontext
+--editable git+https://github.com/majosm/arraycontext.git@empty-subcontainers#egg=arraycontext
 --editable git+https://github.com/inducer/meshmode.git#egg=meshmode
---editable git+https://github.com/inducer/grudge.git#egg=grudge
+--editable git+https://github.com/majosm/grudge.git@use-actx-reductions#egg=grudge
 --editable git+https://github.com/inducer/pytato.git#egg=pytato
 --editable git+https://github.com/ecisneros8/pyrometheus.git#egg=pyrometheus
 --editable git+https://github.com/illinois-ceesd/logpyle.git#egg=logpyle
diff --git a/test/test_euler.py b/test/test_euler.py
index f93cbe70c..db869215c 100644
--- a/test/test_euler.py
+++ b/test/test_euler.py
@@ -53,8 +53,6 @@
     pytest_generate_tests_for_pyopencl_array_context
     as pytest_generate_tests)
 
-from mirgecom.simutil import max_component_norm
-
 from grudge.shortcuts import make_visualizer
 from mirgecom.inviscid import get_inviscid_timestep
 
@@ -235,9 +233,9 @@ def test_vortex_rhs(actx_factory, order):
             discr, eos=IdealSingleGas(), boundaries=boundaries,
             cv=vortex_soln, time=0.0)
 
-        err_max = max_component_norm(discr, inviscid_rhs, np.inf)
+        err_max = discr.norm(inviscid_rhs, np.inf)
 
-        eoc_rec.add_data_point(1.0 / nel_1d, err_max)
+        eoc_rec.add_data_point(1.0 / nel_1d, actx.to_numpy(err_max))
 
     logger.info(
         f"Error for (dim,order) = ({dim},{order}):\n"
@@ -295,7 +293,7 @@ def test_lump_rhs(actx_factory, dim, order):
         )
         expected_rhs = lump.exact_rhs(discr, cv=lump_soln, time=0)
 
-        err_max = max_component_norm(discr, inviscid_rhs-expected_rhs, np.inf)
+        err_max = actx.to_numpy(discr.norm(inviscid_rhs-expected_rhs, np.inf))
         if err_max > maxxerr:
             maxxerr = err_max
 
@@ -505,7 +503,7 @@ def rhs(t, q):
         maxerr = max(write_soln(cv, False))
     else:
         expected_result = initializer(nodes, time=t)
-        maxerr = max_component_norm(discr, cv-expected_result, np.inf)
+        maxerr = actx.to_numpy(discr.norm(cv-expected_result, np.inf))
 
     logger.info(f"Max Error: {maxerr}")
     if maxerr > exittol:
diff --git a/test/test_filter.py b/test/test_filter.py
index 2bbee78f5..f9cb8216b 100644
--- a/test/test_filter.py
+++ b/test/test_filter.py
@@ -30,6 +30,7 @@
 import numpy as np
 from functools import partial
 
+from arraycontext import flatten
 from meshmode.dof_array import thaw
 from grudge.eager import EagerDGDiscretization
 from meshmode.array_context import (  # noqa
@@ -40,6 +41,7 @@
 )
 from meshmode.dof_array import thaw  # noqa
 from mirgecom.filter import make_spectral_filter
+from mirgecom.simutil import componentwise_norms
 
 
 @pytest.mark.parametrize("dim", [1, 2, 3])
@@ -192,8 +194,8 @@ def test_filter_function(actx_factory, dim, order, do_viz=False):
     filtered_soln = filter_modally(discr, "vol", cutoff,
                                    frfunc, uniform_soln)
     soln_resid = uniform_soln - filtered_soln
-    from mirgecom.simutil import componentwise_norms
-    max_errors = componentwise_norms(discr, soln_resid, np.inf)
+    max_errors = actx.to_numpy(
+        flatten(componentwise_norms(discr, soln_resid, np.inf), actx))
 
     tol = 1e-14
 
@@ -218,7 +220,8 @@ def polyfn(coeff):  # , x_vec):
     filtered_field = filter_modally(discr, "vol", cutoff,
                                     frfunc, field)
     soln_resid = field - filtered_field
-    max_errors = [actx.to_numpy(discr.norm(v, np.inf)) for v in soln_resid]
+    max_errors = actx.to_numpy(
+        flatten(componentwise_norms(discr, soln_resid, np.inf), actx))
     logger.info(f"Field = {field}")
     logger.info(f"Filtered = {filtered_field}")
     logger.info(f"Max Errors (poly) = {max_errors}")
@@ -254,6 +257,7 @@ def polyfn(coeff):  # , x_vec):
             ]
             vis.write_vtk_file(f"filter_test_{field_order}.vtu", io_fields)
         field_resid = unfiltered_spectrum - filtered_spectrum
-        max_errors = [actx.to_numpy(discr.norm(v, np.inf)) for v in field_resid]
+        max_errors = actx.to_numpy(
+            flatten(componentwise_norms(discr, field_resid, np.inf), actx))
         # fields should be different, but not too different
         assert(tol > np.max(max_errors) > threshold)
diff --git a/test/test_lazy.py b/test/test_lazy.py
index 5b3472bcc..58da7b71b 100644
--- a/test/test_lazy.py
+++ b/test/test_lazy.py
@@ -76,18 +76,16 @@ def get_discr(order):
 # Mimics math.isclose for state arrays
 def _isclose(discr, x, y, rel_tol=1e-9, abs_tol=0, return_operands=False):
 
-    from mirgecom.simutil import componentwise_norms
-    from arraycontext import flatten
     actx = x.array_context
-    lhs = actx.to_numpy(flatten(componentwise_norms(discr, x - y, np.inf), actx))
+    lhs = discr.norm(x - y, np.inf)
 
-    rhs = np.maximum(
-        rel_tol * np.maximum(
-            actx.to_numpy(flatten(componentwise_norms(discr, x, np.inf), actx)),
-            actx.to_numpy(flatten(componentwise_norms(discr, y, np.inf), actx))),
+    rhs = actx.np.maximum(
+        rel_tol * actx.np.maximum(
+            discr.norm(x, np.inf),
+            discr.norm(y, np.inf)),
         abs_tol)
 
-    is_close = np.all(lhs <= rhs)
+    is_close = actx.to_numpy(actx.np.all(lhs <= rhs))
 
     if return_operands:
         return is_close, lhs, rhs
diff --git a/test/test_operators.py b/test/test_operators.py
index 761d32f97..b8dc0d4ce 100644
--- a/test/test_operators.py
+++ b/test/test_operators.py
@@ -217,14 +217,10 @@ def sym_eval(expr, x_vec):
         test_data = test_func(nodes)
         exact_grad = grad_test_func(nodes)
 
-        from mirgecom.simutil import componentwise_norms
-        from arraycontext import flatten
-
-        err_scale = max(flatten(componentwise_norms(discr, exact_grad, np.inf),
-                                actx))
+        err_scale = discr.norm(exact_grad, np.inf)
 
         if err_scale <= 1e-16:
-            err_scale = 1
+            err_scale = actx.from_numpy(np.array(1))
 
         print(f"{test_data=}")
         print(f"{exact_grad=}")
@@ -236,9 +232,7 @@ def sym_eval(expr, x_vec):
 
         from mirgecom.operators import grad_operator
         test_grad = grad_operator(discr, test_data, test_data_flux_bnd)
-        grad_err = \
-            max(flatten(componentwise_norms(discr, test_grad - exact_grad, np.inf),
-                        actx)) / err_scale
+        grad_err = discr.norm(test_grad - exact_grad, np.inf) / err_scale
 
         eoc.add_data_point(actx.to_numpy(h_max), actx.to_numpy(grad_err))
 
diff --git a/test/test_restart.py b/test/test_restart.py
index 55894c538..71ec62a8d 100644
--- a/test/test_restart.py
+++ b/test/test_restart.py
@@ -76,5 +76,4 @@ def test_restart_cv(actx_factory, nspecies):
     restart_data = read_restart_data(actx, rst_filename)
 
     resid = test_state - restart_data["state"]
-    from mirgecom.simutil import max_component_norm
-    assert max_component_norm(discr, resid, np.inf) == 0
+    assert actx.to_numpy(discr.norm(resid, np.inf)) == 0
diff --git a/test/test_simutil.py b/test/test_simutil.py
index a80cb7fdb..fb8bd0a5c 100644
--- a/test/test_simutil.py
+++ b/test/test_simutil.py
@@ -106,10 +106,9 @@ def test_basic_cfd_healthcheck(actx_factory):
                                  max_value=np.inf)
 
 
-def test_analytic_comparison(actx_factory):
+def test_componentwise_norms(actx_factory):
     """Quick test of state comparison routine."""
-    from mirgecom.initializers import Vortex2D
-    from mirgecom.simutil import compare_fluid_solutions, componentwise_norms
+    from mirgecom.simutil import componentwise_norms
 
     actx = actx_factory()
     nel_1d = 4
@@ -130,17 +129,17 @@ def test_analytic_comparison(actx_factory):
     energy = ones
     velocity = 2 * nodes
     mom = mass * velocity
-    vortex_init = Vortex2D()
-    vortex_soln = vortex_init(x_vec=nodes, eos=IdealSingleGas())
 
     cv = make_conserved(dim, mass=mass, energy=energy, momentum=mom)
-    resid = vortex_soln - cv
 
-    expected_errors = actx.to_numpy(
-        flatten(componentwise_norms(discr, resid, order=np.inf), actx)).tolist()
-
-    errors = compare_fluid_solutions(discr, cv, cv)
-    assert max(errors) == 0
-
-    errors = compare_fluid_solutions(discr, cv, vortex_soln)
-    assert errors == expected_errors
+    component_norms = componentwise_norms(discr, cv)
+
+    assert (
+        actx.to_numpy(component_norms.mass)
+        == actx.to_numpy(discr.norm(cv.mass, np.inf)))
+    assert (
+        actx.to_numpy(component_norms.energy)
+        == actx.to_numpy(discr.norm(cv.energy, np.inf)))
+    assert np.all(
+        actx.to_numpy(flatten(component_norms.momentum, actx))
+        == actx.to_numpy(flatten(discr.norm(cv.momentum, np.inf), actx)))