diff --git a/examples/advection/surface.py b/examples/advection/surface.py
index bc06efc5d5f4324c32d78a2e81cc08847d216180..0e986f31b357c868e115f90a6e982611b22629ad 100644
--- a/examples/advection/surface.py
+++ b/examples/advection/surface.py
@@ -31,7 +31,6 @@ import numpy as np
 import pyopencl as cl
 import pyopencl.tools as cl_tools
 
-from arraycontext import thaw
 from grudge.array_context import PyOpenCLArrayContext
 
 from meshmode.dof_array import flatten
@@ -62,7 +61,7 @@ class Plotter:
             import matplotlib.pyplot as pt
             self.fig = pt.figure(figsize=(8, 8), dpi=300)
 
-            x = thaw(dcoll.discr_from_dd(dof_desc.DD_VOLUME).nodes(), actx)
+            x = actx.thaw(dcoll.discr_from_dd(dof_desc.DD_VOLUME).nodes())
             self.x = actx.to_numpy(flatten(actx.np.arctan2(x[1], x[0])))
         elif self.ambient_dim == 3:
             from grudge.shortcuts import make_visualizer
@@ -174,7 +173,7 @@ def main(ctx_factory, dim=2, order=4, use_quad=False, visualize=False):
     # {{{ Surface advection operator
 
     # velocity field
-    x = thaw(dcoll.nodes(), actx)
+    x = actx.thaw(dcoll.nodes())
     c = make_obj_array([-x[1], x[0], 0.0])[:dim]
 
     def f_initial_condition(x):
@@ -238,7 +237,7 @@ def main(ctx_factory, dim=2, order=4, use_quad=False, visualize=False):
 
         df = dof_desc.DOFDesc(FACE_RESTR_INTERIOR)
         face_discr = dcoll.discr_from_dd(df)
-        face_normal = thaw(dcoll.normal(dd=df), actx)
+        face_normal = actx.thaw(dcoll.normal(dd=df))
 
         from meshmode.discretization.visualization import make_visualizer
         vis = make_visualizer(actx, face_discr)
diff --git a/examples/advection/var-velocity.py b/examples/advection/var-velocity.py
index 289f1274de7b1e6659c23d491cd217f27da535a8..de1b45354467e2c15ad98c0b41be0f7a21c58b20 100644
--- a/examples/advection/var-velocity.py
+++ b/examples/advection/var-velocity.py
@@ -29,7 +29,6 @@ import numpy as np
 import pyopencl as cl
 import pyopencl.tools as cl_tools
 
-from arraycontext import thaw
 from grudge.array_context import PyOpenCLArrayContext
 
 from meshmode.dof_array import flatten
@@ -61,7 +60,7 @@ class Plotter:
             self.ylim = ylim
 
             volume_discr = dcoll.discr_from_dd(dof_desc.DD_VOLUME)
-            self.x = actx.to_numpy(flatten(thaw(volume_discr.nodes()[0], actx)))
+            self.x = actx.to_numpy(flatten(actx.thaw(volume_discr.nodes()[0])))
         else:
             from grudge.shortcuts import make_visualizer
             self.vis = make_visualizer(dcoll)
@@ -168,7 +167,7 @@ def main(ctx_factory, dim=2, order=4, use_quad=False, visualize=False,
 
     from grudge.models.advection import VariableCoefficientAdvectionOperator
 
-    x = thaw(dcoll.nodes(), actx)
+    x = actx.thaw(dcoll.nodes())
 
     # velocity field
     if dim == 1:
diff --git a/examples/advection/weak.py b/examples/advection/weak.py
index 433f9092e9276e7d9fdcc9c5b5079c2088468960..3470fdd607479ff11b04401c76b7e3d148a3c9f1 100644
--- a/examples/advection/weak.py
+++ b/examples/advection/weak.py
@@ -30,7 +30,6 @@ import numpy.linalg as la
 import pyopencl as cl
 import pyopencl.tools as cl_tools
 
-from arraycontext import thaw
 from grudge.array_context import PyOpenCLArrayContext
 
 from meshmode.dof_array import flatten
@@ -60,7 +59,7 @@ class Plotter:
             self.ylim = ylim
 
             volume_discr = dcoll.discr_from_dd(dof_desc.DD_VOLUME)
-            self.x = actx.to_numpy(flatten(thaw(volume_discr.nodes()[0], actx)))
+            self.x = actx.to_numpy(flatten(actx.thaw(volume_discr.nodes()[0])))
         else:
             from grudge.shortcuts import make_visualizer
             self.vis = make_visualizer(dcoll)
@@ -152,13 +151,13 @@ def main(ctx_factory, dim=2, order=4, visualize=False):
         dcoll,
         c,
         inflow_u=lambda t: u_analytic(
-            thaw(dcoll.nodes(dd=BTAG_ALL), actx),
+            actx.thaw(dcoll.nodes(dd=BTAG_ALL)),
             t=t
         ),
         flux_type=flux_type
     )
 
-    nodes = thaw(dcoll.nodes(), actx)
+    nodes = actx.thaw(dcoll.nodes())
     u = u_analytic(nodes, t=0)
 
     def rhs(t, u):
diff --git a/examples/euler/acoustic_pulse.py b/examples/euler/acoustic_pulse.py
index 2672dda127eaec7b8bd8143f75e70d2d4de118a2..779062910558c8e011353c9a7fc0422e5b84fa20 100644
--- a/examples/euler/acoustic_pulse.py
+++ b/examples/euler/acoustic_pulse.py
@@ -28,7 +28,6 @@ import numpy as np
 import pyopencl as cl
 import pyopencl.tools as cl_tools
 
-from arraycontext import thaw, freeze
 from grudge.array_context import (
     PyOpenCLArrayContext,
     PytatoPyOpenCLArrayContext
@@ -175,7 +174,7 @@ def run_acoustic_pulse(actx,
     cn = 0.5*(order + 1)**2
     dt = cfl * actx.to_numpy(h_min_from_volume(dcoll)) / cn
 
-    fields = acoustic_pulse_condition(thaw(dcoll.nodes(), actx))
+    fields = acoustic_pulse_condition(actx.thaw(dcoll.nodes()))
 
     logger.info("Timestep size: %g", dt)
 
@@ -204,7 +203,7 @@ def run_acoustic_pulse(actx,
                 )
             assert norm_q < 5
 
-        fields = thaw(freeze(fields, actx), actx)
+        fields = actx.thaw(actx.freeze(fields))
         fields = rk4_step(fields, t, dt, compiled_rhs)
         t += dt
         step += 1
diff --git a/examples/euler/vortex.py b/examples/euler/vortex.py
index 0e4cf7d7691baadbb09fac53474bf19f4a025935..9f00743e53747c045801b7dda6e9da89fbaf38ad 100644
--- a/examples/euler/vortex.py
+++ b/examples/euler/vortex.py
@@ -26,8 +26,6 @@ THE SOFTWARE.
 import pyopencl as cl
 import pyopencl.tools as cl_tools
 
-from arraycontext import thaw, freeze
-
 from grudge.array_context import PytatoPyOpenCLArrayContext, PyOpenCLArrayContext
 from grudge.models.euler import (
     vortex_initial_condition,
@@ -111,7 +109,7 @@ def run_vortex(actx, order=3, resolution=8, final_time=5,
 
     compiled_rhs = actx.compile(rhs)
 
-    fields = vortex_initial_condition(thaw(dcoll.nodes(), actx))
+    fields = vortex_initial_condition(actx.thaw(dcoll.nodes()))
 
     from grudge.dt_utils import h_min_from_volume
 
@@ -147,7 +145,7 @@ def run_vortex(actx, order=3, resolution=8, final_time=5,
                 )
             assert norm_q < 200
 
-        fields = thaw(freeze(fields, actx), actx)
+        fields = actx.thaw(actx.freeze(fields))
         fields = rk4_step(fields, t, dt, compiled_rhs)
         t += dt
         step += 1
diff --git a/examples/geometry.py b/examples/geometry.py
index f79cd6bb904e7c40599bc70921369f10de8295fe..442bbcfff4a94d5b1ffe59928c1544de591e078e 100644
--- a/examples/geometry.py
+++ b/examples/geometry.py
@@ -30,7 +30,6 @@ import numpy as np  # noqa
 import pyopencl as cl
 import pyopencl.tools as cl_tools
 
-from arraycontext import thaw
 from grudge.array_context import PyOpenCLArrayContext
 
 from grudge import DiscretizationCollection, shortcuts
@@ -51,9 +50,9 @@ def main(write_output=True):
 
     dcoll = DiscretizationCollection(actx, mesh, order=4)
 
-    nodes = thaw(dcoll.nodes(), actx)
-    bdry_nodes = thaw(dcoll.nodes(dd=BTAG_ALL), actx)
-    bdry_normals = thaw(dcoll.normal(dd=BTAG_ALL), actx)
+    nodes = actx.thaw(dcoll.nodes())
+    bdry_nodes = actx.thaw(dcoll.nodes(dd=BTAG_ALL))
+    bdry_normals = actx.thaw(dcoll.normal(dd=BTAG_ALL))
 
     if write_output:
         vis = shortcuts.make_visualizer(dcoll)
diff --git a/examples/hello-grudge.py b/examples/hello-grudge.py
index bf93fd08b93eab03d854f6e14652e8ae1fa60e83..cfb724115a487ec8c583184d3ee06ad2af4d1121 100644
--- a/examples/hello-grudge.py
+++ b/examples/hello-grudge.py
@@ -14,7 +14,6 @@ from grudge.discretization import DiscretizationCollection
 import grudge.op as op
 from meshmode.mesh.generation import generate_box_mesh
 from meshmode.array_context import PyOpenCLArrayContext
-from arraycontext import thaw
 from grudge.dof_desc import DTAG_BOUNDARY, FACE_RESTR_INTERIOR
 
 
@@ -43,7 +42,7 @@ def left_boundary_condition(x, t):
 def flux(dcoll, u_tpair):
     dd = u_tpair.dd
     velocity = np.array([2 * np.pi])
-    normal = thaw(dcoll.normal(dd), actx)
+    normal = actx.thaw(dcoll.normal(dd))
 
     v_dot_n = np.dot(velocity, normal)
     u_upwind = actx.np.where(v_dot_n > 0,
@@ -55,8 +54,8 @@ vol_discr = dcoll.discr_from_dd("vol")
 left_bndry = DTAG_BOUNDARY("left")
 right_bndry = DTAG_BOUNDARY("right")
 
-x_vol = thaw(dcoll.nodes(), actx)
-x_bndry = thaw(dcoll.discr_from_dd(left_bndry).nodes(), actx)
+x_vol = actx.thaw(dcoll.nodes())
+x_bndry = actx.thaw(dcoll.discr_from_dd(left_bndry).nodes())
 
 uh = initial_condition(x_vol)
 
diff --git a/examples/maxwell/cavities.py b/examples/maxwell/cavities.py
index 90a515bdffaedec4e5ceca124532ca391d0a29c3..3d581c18a9d20378977e9b85b73b58d795ec7667 100644
--- a/examples/maxwell/cavities.py
+++ b/examples/maxwell/cavities.py
@@ -28,7 +28,6 @@ import numpy as np
 import pyopencl as cl
 import pyopencl.tools as cl_tools
 
-from arraycontext import thaw
 from grudge.array_context import PyOpenCLArrayContext
 
 from grudge.shortcuts import set_up_rk4
@@ -84,7 +83,7 @@ def main(ctx_factory, dim=3, order=4, visualize=False):
         else:
             return get_rectangular_cavity_mode(actx, x, t, 1, (2, 3))
 
-    fields = cavity_mode(thaw(dcoll.nodes(), actx), t=0)
+    fields = cavity_mode(actx.thaw(dcoll.nodes()), t=0)
 
     maxwell_operator.check_bc_coverage(mesh)
 
diff --git a/examples/wave/var-propagation-speed.py b/examples/wave/var-propagation-speed.py
index 407dd06b266dbe868dac187c37bec2837b67de37..9929f6dbfd6a387da20162b81c9b046682abbf06 100644
--- a/examples/wave/var-propagation-speed.py
+++ b/examples/wave/var-propagation-speed.py
@@ -28,7 +28,6 @@ import numpy as np
 import pyopencl as cl
 import pyopencl.tools as cl_tools
 
-from arraycontext import thaw
 from grudge.array_context import PyOpenCLArrayContext
 
 from grudge.shortcuts import set_up_rk4
@@ -63,7 +62,7 @@ def main(ctx_factory, dim=2, order=4, visualize=False):
         source_center = np.array([0.1, 0.22, 0.33])[:dcoll.dim]
         source_width = 0.05
         source_omega = 3
-        nodes = thaw(dcoll.nodes(), actx)
+        nodes = actx.thaw(dcoll.nodes())
         source_center_dist = flat_obj_array(
             [nodes[i] - source_center[i] for i in range(dcoll.dim)]
         )
@@ -75,7 +74,7 @@ def main(ctx_factory, dim=2, order=4, visualize=False):
             )
         )
 
-    x = thaw(dcoll.nodes(), actx)
+    x = actx.thaw(dcoll.nodes())
     ones = dcoll.zeros(actx) + 1
     c = actx.np.where(np.dot(x, x) < 0.15, 0.1 * ones, 0.2 * ones)
 
diff --git a/examples/wave/wave-min-mpi.py b/examples/wave/wave-min-mpi.py
index c58ed841d30b4d8da81e4516fe3b92ef2ba69093..6c56353bdcac4fb044433b8148646dc203e825d2 100644
--- a/examples/wave/wave-min-mpi.py
+++ b/examples/wave/wave-min-mpi.py
@@ -28,7 +28,6 @@ import numpy as np
 import pyopencl as cl
 import pyopencl.tools as cl_tools
 
-from arraycontext import thaw
 from grudge.array_context import MPIPyOpenCLArrayContext
 
 from grudge.shortcuts import set_up_rk4
@@ -88,7 +87,7 @@ def main(ctx_factory, dim=2, order=4, visualize=False):
         source_center = np.array([0.1, 0.22, 0.33])[:dcoll.dim]
         source_width = 0.05
         source_omega = 3
-        nodes = thaw(dcoll.nodes(), actx)
+        nodes = actx.thaw(dcoll.nodes())
         source_center_dist = flat_obj_array(
             [nodes[i] - source_center[i] for i in range(dcoll.dim)]
         )
diff --git a/examples/wave/wave-op-mpi.py b/examples/wave/wave-op-mpi.py
index e332404436d96544b6aaf529a51573ee1feac860..8c23336d03a3c226c4284da05ebad31328d7f578 100644
--- a/examples/wave/wave-op-mpi.py
+++ b/examples/wave/wave-op-mpi.py
@@ -30,7 +30,6 @@ import pyopencl as cl
 import pyopencl.tools as cl_tools
 
 from arraycontext import (
-    thaw,
     with_container_arithmetic,
     dataclass_array_container
 )
@@ -73,12 +72,12 @@ class WaveState:
         return self.u.array_context
 
 
-def wave_flux(dcoll, c, w_tpair):
+def wave_flux(actx, dcoll, c, w_tpair):
     u = w_tpair.u
     v = w_tpair.v
     dd = w_tpair.dd
 
-    normal = thaw(dcoll.normal(dd), u.int.array_context)
+    normal = actx.thaw(dcoll.normal(dd))
 
     flux_weak = WaveState(
         u=v.avg @ normal,
@@ -99,7 +98,7 @@ class _WaveStateTag:
     pass
 
 
-def wave_operator(dcoll, c, w, quad_tag=None):
+def wave_operator(actx, dcoll, c, w, quad_tag=None):
     dd_base = as_dofdesc("vol")
     dd_vol = DOFDesc("vol", quad_tag)
     dd_faces = DOFDesc("all_faces", quad_tag)
@@ -135,13 +134,14 @@ def wave_operator(dcoll, c, w, quad_tag=None):
                 dcoll,
                 dd_faces,
                 wave_flux(
+                    actx,
                     dcoll, c=c,
                     w_tpair=op.bdry_trace_pair(dcoll,
                                                dd_btag,
                                                interior=dir_bval,
                                                exterior=dir_bc)
                 ) + sum(
-                    wave_flux(dcoll, c=c, w_tpair=interp_to_surf_quad(tpair))
+                    wave_flux(actx, dcoll, c=c, w_tpair=interp_to_surf_quad(tpair))
                     for tpair in op.interior_trace_pairs(dcoll, w,
                         comm_tag=_WaveStateTag)
                 )
@@ -165,7 +165,7 @@ def bump(actx, dcoll, t=0):
     source_width = 0.05
     source_omega = 3
 
-    nodes = thaw(dcoll.nodes(), actx)
+    nodes = actx.thaw(dcoll.nodes())
     center_dist = flat_obj_array([
         nodes[i] - source_center[i]
         for i in range(dcoll.dim)
@@ -258,7 +258,7 @@ def main(ctx_factory, dim=2, order=3,
     vis = make_visualizer(dcoll)
 
     def rhs(t, w):
-        return wave_operator(dcoll, c=c, w=w, quad_tag=quad_tag)
+        return wave_operator(actx, dcoll, c=c, w=w, quad_tag=quad_tag)
 
     compiled_rhs = actx.compile(rhs)
 
diff --git a/examples/wave/wave-op-var-velocity.py b/examples/wave/wave-op-var-velocity.py
index 9dace53c53d1e1c11bc8424bd34a298493a73bcc..43c72eff9ff2242c3265386e88512acf50a43517 100644
--- a/examples/wave/wave-op-var-velocity.py
+++ b/examples/wave/wave-op-var-velocity.py
@@ -29,7 +29,6 @@ import numpy.linalg as la  # noqa
 import pyopencl as cl
 import pyopencl.tools as cl_tools
 
-from arraycontext import thaw
 from grudge.array_context import PyOpenCLArrayContext
 
 from pytools.obj_array import flat_obj_array
@@ -48,14 +47,14 @@ logger = logging.getLogger(__name__)
 
 # {{{ wave equation bits
 
-def wave_flux(dcoll, c, w_tpair):
+def wave_flux(actx, dcoll, c, w_tpair):
     dd = w_tpair.dd
     dd_quad = dd.with_discr_tag(DISCR_TAG_QUAD)
 
     u = w_tpair[0]
     v = w_tpair[1:]
 
-    normal = thaw(dcoll.normal(dd), u.int.array_context)
+    normal = actx.thaw(dcoll.normal(dd))
 
     flux_weak = flat_obj_array(
             np.dot(v.avg, normal),
@@ -76,7 +75,7 @@ def wave_flux(dcoll, c, w_tpair):
     return op.project(dcoll, dd_quad, dd_allfaces_quad, c_quad*flux_quad)
 
 
-def wave_operator(dcoll, c, w):
+def wave_operator(actx, dcoll, c, w):
     u = w[0]
     v = w[1:]
 
@@ -104,13 +103,14 @@ def wave_operator(dcoll, c, w):
                 dcoll,
                 dd_allfaces_quad,
                 wave_flux(
+                    actx,
                     dcoll, c=c,
                     w_tpair=op.bdry_trace_pair(dcoll,
                                                BTAG_ALL,
                                                interior=dir_bval,
                                                exterior=dir_bc)
                 ) + sum(
-                    wave_flux(dcoll, c=c, w_tpair=tpair)
+                    wave_flux(actx, dcoll, c=c, w_tpair=tpair)
                     for tpair in op.interior_trace_pairs(dcoll, w)
                 )
             )
@@ -135,7 +135,7 @@ def bump(actx, dcoll, t=0, width=0.05, center=None):
     center = center[:dcoll.dim]
     source_omega = 3
 
-    nodes = thaw(dcoll.nodes(), actx)
+    nodes = actx.thaw(dcoll.nodes())
     center_dist = flat_obj_array([
         nodes[i] - center[i]
         for i in range(dcoll.dim)
@@ -189,7 +189,7 @@ def main(ctx_factory, dim=2, order=3, visualize=False):
     vis = make_visualizer(dcoll)
 
     def rhs(t, w):
-        return wave_operator(dcoll, c=c, w=w)
+        return wave_operator(actx, dcoll, c=c, w=w)
 
     logger.info("dt = %g", dt)
 
diff --git a/grudge/array_context.py b/grudge/array_context.py
index f1a27ba3e8b32b6a7d72dd4bcefb3339bfcdcd06..fb6788359147b8a65264af68af1129ec9ee9d438 100644
--- a/grudge/array_context.py
+++ b/grudge/array_context.py
@@ -176,7 +176,7 @@ class _DistributedLazilyCompilingFunctionCaller(
                 part_id_to_prg[part.pid],
                 part_prg_name_to_tags,
                 part_prg_name_to_axes
-            ) = self._dag_to_transformed_loopy_prg(d)
+            ) = self._dag_to_transformed_pytato_prg(d)
 
             assert not (set(name_in_program_to_tags.keys())
                         & set(part_prg_name_to_tags.keys()))
diff --git a/grudge/discretization.py b/grudge/discretization.py
index 3a86b70603a7f018f62dd6a238e40f32a1f6657a..43bd242266688091519a5087fa21f600b4ab09df 100644
--- a/grudge/discretization.py
+++ b/grudge/discretization.py
@@ -734,10 +734,9 @@ class DiscretizationCollection:
         :arg dd: a :class:`~grudge.dof_desc.DOFDesc` as the surface discretization.
         :returns: an object array of frozen :class:`~meshmode.dof_array.DOFArray`\ s.
         """
-        from arraycontext import freeze
         from grudge.geometry import normal
 
-        return freeze(normal(self._setup_actx, self, dd))
+        return self._setup_actx.freeze(normal(self._setup_actx, self, dd))
 
     # }}}
 
diff --git a/grudge/dt_utils.py b/grudge/dt_utils.py
index a9fd8cf43ae3c5b3239b51658fa6093dfa5aa5a3..afda657f01f0994e2b76d465dc4a43256aa94e2c 100644
--- a/grudge/dt_utils.py
+++ b/grudge/dt_utils.py
@@ -45,7 +45,7 @@ THE SOFTWARE.
 
 import numpy as np
 
-from arraycontext import ArrayContext, thaw, freeze, Scalar, tag_axes
+from arraycontext import ArrayContext, Scalar, tag_axes
 from meshmode.transform_metadata import (FirstAxisIsElementsTag,
                                          DiscretizationDOFAxisTag,
                                          DiscretizationFaceAxisTag,
@@ -93,7 +93,7 @@ def characteristic_lengthscales(
     @memoize_in(dcoll, (characteristic_lengthscales,
                         "compute_characteristic_lengthscales"))
     def _compute_characteristic_lengthscales():
-        return freeze(
+        return actx.freeze(
             DOFArray(
                 actx,
                 data=tuple(
@@ -102,12 +102,12 @@ def characteristic_lengthscales(
                     cng * geo_facts
                     for cng, geo_facts in zip(
                         dt_non_geometric_factors(dcoll),
-                        thaw(dt_geometric_factors(dcoll), actx)
+                        actx.thaw(dt_geometric_factors(dcoll))
                     )
                 )
             )
         )
-    return thaw(_compute_characteristic_lengthscales(), actx)
+    return actx.thaw(_compute_characteristic_lengthscales())
 
 
 @memoize_on_first_arg
@@ -270,7 +270,7 @@ def dt_geometric_factors(
 
     if dcoll.dim == 1:
         # Inscribed "circle" radius is half the cell size
-        return freeze(cell_vols/2)
+        return actx.freeze(cell_vols/2)
 
     dd_face = DOFDesc("all_faces", dd.discretization_tag)
     face_discr = dcoll.discr_from_dd(dd_face)
@@ -327,7 +327,7 @@ def dt_geometric_factors(
             )
         )
 
-    return freeze(DOFArray(
+    return actx.freeze(DOFArray(
         actx,
         data=tuple(
             actx.einsum("e,ei->ei", 1/sae_i, cv_i,
diff --git a/grudge/geometry/metrics.py b/grudge/geometry/metrics.py
index d065a4360e1830149130859868a25becbd29b7b7..9110201be25cd69e11c8fe0eefbdb3c5dcf39ad8 100644
--- a/grudge/geometry/metrics.py
+++ b/grudge/geometry/metrics.py
@@ -60,7 +60,7 @@ THE SOFTWARE.
 
 import numpy as np
 
-from arraycontext import thaw, freeze, ArrayContext, tag_axes
+from arraycontext import ArrayContext, tag_axes
 from meshmode.dof_array import DOFArray
 
 from grudge import DiscretizationCollection
@@ -173,7 +173,7 @@ def forward_metric_nth_derivative(
     vec = num_reference_derivative(
         dcoll.discr_from_dd(inner_dd),
         flat_ref_axes,
-        thaw(dcoll.discr_from_dd(inner_dd).nodes(), actx)[xyz_axis]
+        actx.thaw(dcoll.discr_from_dd(inner_dd).nodes())[xyz_axis]
     )
 
     return _geometry_to_quad_if_requested(
@@ -526,12 +526,12 @@ def inverse_surface_metric_derivative_mat(
                 for rst_axis in range(dcoll.dim)])
             for xyz_axis in range(dcoll.ambient_dim)])
 
-        return freeze(tag_axes(actx, {
+        return actx.freeze(tag_axes(actx, {
             0: DiscretizationAmbientDimAxisTag(),
             1: DiscretizationTopologicalDimAxisTag()},
             mat))
 
-    return thaw(_inv_surf_metric_deriv(), actx)
+    return actx.thaw(_inv_surf_metric_deriv())
 
 
 def _signed_face_ones(
@@ -549,13 +549,13 @@ def _signed_face_ones(
         actx, dtype=dcoll.real_dtype
     ) + 1
 
-    from arraycontext import to_numpy, from_numpy, thaw
+    from arraycontext import to_numpy, from_numpy
 
     _signed_face_ones_numpy = to_numpy(signed_ones, actx)
 
     for igrp, grp in enumerate(all_faces_conn.groups):
         for batch in grp.batches:
-            i = to_numpy(thaw(batch.to_element_indices, actx), actx)
+            i = to_numpy(actx.thaw(batch.to_element_indices), actx)
             grp_field = _signed_face_ones_numpy[igrp].reshape(-1)
             grp_field[i] = \
                 (2.0 * (batch.to_element_face % 2) - 1.0) * grp_field[i]
@@ -644,9 +644,9 @@ def area_element(
                 actx, dcoll, dd=dd,
                 _use_geoderiv_connection=_use_geoderiv_connection).norm_squared())
 
-        return freeze(result, actx)
+        return actx.freeze(result)
 
-    return thaw(_area_elements(), actx)
+    return actx.thaw(_area_elements())
 
 # }}}
 
@@ -751,10 +751,9 @@ def mv_normal(
 
             result = mv / actx.np.sqrt(mv.norm_squared())
 
-        return freeze(result, actx)
+        return actx.freeze(result)
 
-    n = _normal()
-    return thaw(n, actx)
+    return actx.thaw(_normal())
 
 
 def normal(actx: ArrayContext, dcoll: DiscretizationCollection, dd,
diff --git a/grudge/models/advection.py b/grudge/models/advection.py
index aacd2fcfdf1b80f119900cf0e22f0fe4e2ef21b3..cfe1a492015c02adc02b122c916ae82f4d1bc607 100644
--- a/grudge/models/advection.py
+++ b/grudge/models/advection.py
@@ -30,8 +30,6 @@ import numpy as np
 import grudge.op as op
 import types
 
-from arraycontext.container.traversal import thaw
-
 from grudge.models import HyperbolicOperator
 
 
@@ -43,7 +41,7 @@ def advection_weak_flux(dcoll, flux_type, u_tpair, velocity):
     """
     actx = u_tpair.int.array_context
     dd = u_tpair.dd
-    normal = thaw(dcoll.normal(dd), actx)
+    normal = actx.thaw(dcoll.normal(dd))
     v_dot_n = np.dot(velocity, normal)
 
     flux_type = flux_type.lower()
@@ -92,7 +90,7 @@ class StrongAdvectionOperator(AdvectionOperatorBase):
     def flux(self, u_tpair):
         actx = u_tpair.int.array_context
         dd = u_tpair.dd
-        normal = thaw(self.dcoll.normal(dd), actx)
+        normal = actx.thaw(self.dcoll.normal(dd))
         v_dot_normal = np.dot(self.v, normal)
 
         return u_tpair.int * v_dot_normal - self.weak_flux(u_tpair)
@@ -285,7 +283,7 @@ def v_dot_n_tpair(actx, dcoll, velocity, trace_dd):
     from grudge.trace_pair import TracePair
     from meshmode.discretization.connection import FACE_RESTR_INTERIOR
 
-    normal = thaw(dcoll.normal(trace_dd.with_discr_tag(None)), actx)
+    normal = actx.thaw(dcoll.normal(trace_dd.with_discr_tag(None)))
     v_dot_n = velocity.dot(normal)
     i = op.project(dcoll, trace_dd.with_discr_tag(None), trace_dd, v_dot_n)
 
diff --git a/grudge/models/em.py b/grudge/models/em.py
index e14341cfa5f9bbe3635e5fa1476a935538611df4..7bc9524379e70532ce1ede2b832f9768c4ed813f 100644
--- a/grudge/models/em.py
+++ b/grudge/models/em.py
@@ -28,7 +28,7 @@ THE SOFTWARE.
 """
 
 
-from arraycontext import thaw, get_container_context_recursively
+from arraycontext import get_container_context_recursively
 
 from grudge.models import HyperbolicOperator
 
@@ -121,7 +121,7 @@ class MaxwellOperator(HyperbolicOperator):
         """
 
         actx = get_container_context_recursively(wtpair)
-        normal = thaw(self.dcoll.normal(wtpair.dd), actx)
+        normal = actx.thaw(self.dcoll.normal(wtpair.dd))
 
         if self.fixed_material:
             e, h = self.split_eh(wtpair)
@@ -222,7 +222,7 @@ class MaxwellOperator(HyperbolicOperator):
         """
 
         actx = get_container_context_recursively(w)
-        absorb_normal = thaw(self.dcoll.normal(dd=self.absorb_tag), actx)
+        absorb_normal = actx.thaw(self.dcoll.normal(dd=self.absorb_tag))
 
         e, h = self.split_eh(w)
 
diff --git a/grudge/models/euler.py b/grudge/models/euler.py
index 98acbc1ef5852792b0adcee3fb8d2b77d5ad98ff..f4d6f8f4c98410131459b22817a83374dc1b2209 100644
--- a/grudge/models/euler.py
+++ b/grudge/models/euler.py
@@ -50,7 +50,6 @@ from abc import ABCMeta, abstractmethod
 
 from dataclasses import dataclass
 from arraycontext import (
-    thaw,
     dataclass_array_container,
     with_container_arithmetic
 )
@@ -191,7 +190,7 @@ class PrescribedBC(InviscidBCObject):
         return TracePair(
             dd_bc,
             interior=op.project(dcoll, dd_base, dd_bc, state),
-            exterior=self.prescribed_state(thaw(dcoll.nodes(dd_bc), actx), t=t)
+            exterior=self.prescribed_state(actx.thaw(dcoll.nodes(dd_bc)), t=t)
         )
 
 
@@ -204,7 +203,7 @@ class InviscidWallBC(InviscidBCObject):
             state: ConservedEulerField, t=0):
         actx = state.array_context
         dd_base = as_dofdesc("vol").with_discr_tag(DISCR_TAG_BASE)
-        nhat = thaw(dcoll.normal(dd_bc), actx)
+        nhat = actx.thaw(dcoll.normal(dd_bc))
         interior = op.project(dcoll, dd_base, dd_bc, state)
 
         return TracePair(
@@ -271,7 +270,7 @@ def euler_numerical_flux(
         exterior=euler_volume_flux(dcoll, q_rr, gamma=gamma)
     )
     num_flux = flux_tpair.avg
-    normal = thaw(dcoll.normal(dd_intfaces), actx)
+    normal = actx.thaw(dcoll.normal(dd_intfaces))
 
     if lf_stabilization:
         from arraycontext import outer
diff --git a/grudge/models/wave.py b/grudge/models/wave.py
index 209575ba3258f5cfcf983dc49b13f880ef2c8748..ff8ee57e8151717b7628833b0ecb0ffe9394bbf2 100644
--- a/grudge/models/wave.py
+++ b/grudge/models/wave.py
@@ -28,8 +28,6 @@ THE SOFTWARE.
 
 import numpy as np
 
-from arraycontext import thaw, freeze
-
 from grudge.models import HyperbolicOperator
 
 from meshmode.mesh import BTAG_ALL, BTAG_NONE
@@ -93,7 +91,7 @@ class WeakWaveOperator(HyperbolicOperator):
         u = wtpair[0]
         v = wtpair[1:]
         actx = u.int.array_context
-        normal = thaw(self.dcoll.normal(wtpair.dd), actx)
+        normal = actx.thaw(self.dcoll.normal(wtpair.dd))
 
         central_flux_weak = -self.c*flat_obj_array(
                 np.dot(v.avg, normal),
@@ -136,7 +134,7 @@ class WeakWaveOperator(HyperbolicOperator):
         neu_bc = flat_obj_array(neu_u, -neu_v)
 
         # radiation BCs -------------------------------------------------------
-        rad_normal = thaw(dcoll.normal(dd=self.radiation_tag), actx)
+        rad_normal = actx.thaw(dcoll.normal(dd=self.radiation_tag))
 
         rad_u = op.project(dcoll, "vol", self.radiation_tag, u)
         rad_v = op.project(dcoll, "vol", self.radiation_tag, v)
@@ -214,7 +212,7 @@ class VariableCoefficientWeakWaveOperator(HyperbolicOperator):
             neumann_tag=BTAG_NONE,
             radiation_tag=BTAG_NONE):
         """
-        :arg c: a thawed (with *actx*) :class:`~meshmode.dof_array.DOFArray`
+        :arg c: a frozen :class:`~meshmode.dof_array.DOFArray`
             representing the propogation speed of the wave.
         """
 
@@ -223,11 +221,13 @@ class VariableCoefficientWeakWaveOperator(HyperbolicOperator):
 
         actx = c.array_context
         self.dcoll = dcoll
-        self.c = freeze(c)
+        self.c = c
         self.source_f = source_f
 
         ones = dcoll.zeros(actx) + 1
-        self.sign = freeze(actx.np.where(actx.np.greater(c, 0), ones, -ones))
+        thawed_c = dcoll._setup_actx.thaw(c)
+        self.sign = dcoll._setup_actx.freeze(
+                actx.np.where(actx.np.greater(thawed_c, 0), ones, -ones))
 
         self.dirichlet_tag = dirichlet_tag
         self.neumann_tag = neumann_tag
@@ -242,7 +242,7 @@ class VariableCoefficientWeakWaveOperator(HyperbolicOperator):
         u = wtpair[1]
         v = wtpair[2:]
         actx = u.int.array_context
-        normal = thaw(self.dcoll.normal(wtpair.dd), actx)
+        normal = actx.thaw(self.dcoll.normal(wtpair.dd))
 
         flux_central_weak = -0.5 * flat_obj_array(
             np.dot(v.int*c.int + v.ext*c.ext, normal),
@@ -266,7 +266,7 @@ class VariableCoefficientWeakWaveOperator(HyperbolicOperator):
         v = w[1:]
         actx = u.array_context
 
-        c = thaw(self.c, actx)
+        c = actx.thaw(self.c)
 
         flux_w = flat_obj_array(c, w)
 
@@ -294,13 +294,12 @@ class VariableCoefficientWeakWaveOperator(HyperbolicOperator):
         neu_bc = flat_obj_array(neu_c, neu_u, -neu_v)
 
         # radiation BCs -------------------------------------------------------
-        rad_normal = thaw(dcoll.normal(dd=self.radiation_tag), actx)
+        rad_normal = actx.thaw(dcoll.normal(dd=self.radiation_tag))
 
         rad_c = op.project(dcoll, "vol", self.radiation_tag, c)
         rad_u = op.project(dcoll, "vol", self.radiation_tag, u)
         rad_v = op.project(dcoll, "vol", self.radiation_tag, v)
-        rad_sign = op.project(dcoll, "vol", self.radiation_tag,
-                thaw(self.sign, actx))
+        rad_sign = op.project(dcoll, "vol", self.radiation_tag, actx.thaw(self.sign))
 
         rad_bc = flat_obj_array(
             rad_c,
@@ -345,7 +344,7 @@ class VariableCoefficientWeakWaveOperator(HyperbolicOperator):
             self.radiation_tag])
 
     def max_characteristic_velocity(self, actx, **kwargs):
-        return actx.np.abs(thaw(self.c, actx))
+        return actx.np.abs(actx.thaw(self.c))
 
 # }}}
 
diff --git a/test/test_dt_utils.py b/test/test_dt_utils.py
index 63eb4b4c3643112e182ad458fc97e0c8ff46db32..95d7e4e350018cd1167936fa5bb04455688f387e 100644
--- a/test/test_dt_utils.py
+++ b/test/test_dt_utils.py
@@ -24,8 +24,6 @@ THE SOFTWARE.
 
 import numpy as np
 
-from arraycontext import thaw
-
 from grudge.array_context import (
     PytestPyOpenCLArrayContextFactory,
     PytestPytatoPyOpenCLArrayContextFactory
@@ -76,7 +74,7 @@ def test_geometric_factors_regular_refinement(actx_factory, name):
         dcoll = DiscretizationCollection(actx, mesh, order=builder.order)
         min_factors.append(
             actx.to_numpy(
-                op.nodal_min(dcoll, "vol", thaw(dt_geometric_factors(dcoll), actx)))
+                op.nodal_min(dcoll, "vol", actx.thaw(dt_geometric_factors(dcoll))))
         )
 
     # Resolution is doubled each refinement, so the ratio of consecutive
@@ -88,7 +86,7 @@ def test_geometric_factors_regular_refinement(actx_factory, name):
     # Make sure it works with empty meshes
     mesh = builder.get_mesh(0, builder.mesh_order)
     dcoll = DiscretizationCollection(actx, mesh, order=builder.order)
-    factors = thaw(dt_geometric_factors(dcoll), actx)  # noqa: F841
+    factors = actx.thaw(dt_geometric_factors(dcoll))  # noqa: F841
 
 
 @pytest.mark.parametrize("name", ["interval", "box2d", "box3d"])
diff --git a/test/test_euler_model.py b/test/test_euler_model.py
index 3f5174e261d3b8905507a9b054e63019ac5c142f..13a479d48f448a0c307cf1958cbd11c224238078 100644
--- a/test/test_euler_model.py
+++ b/test/test_euler_model.py
@@ -27,7 +27,6 @@ import pytest
 from grudge.array_context import PytestPyOpenCLArrayContextFactory
 from arraycontext import (
     pytest_generate_tests_for_array_contexts,
-    thaw, freeze
 )
 pytest_generate_tests = pytest_generate_tests_for_array_contexts(
         [PytestPyOpenCLArrayContextFactory])
@@ -82,7 +81,7 @@ def test_euler_vortex_convergence(actx_factory, order):
             discr_tag_to_group_factory=discr_tag_to_group_factory
         )
         h_max = actx.to_numpy(h_max_from_volume(dcoll, dim=dcoll.ambient_dim))
-        nodes = thaw(dcoll.nodes(), actx)
+        nodes = actx.thaw(dcoll.nodes())
 
         # }}}
 
@@ -115,7 +114,7 @@ def test_euler_vortex_convergence(actx_factory, order):
         t = 0.0
         last_q = None
         while t < final_time:
-            fields = thaw(freeze(fields, actx), actx)
+            fields = actx.thaw(actx.freeze(fields))
             fields = rk4_step(fields, t, dt, compiled_rhs)
             t += dt
             logger.info("[%04d] t = %.5f", step, t)
diff --git a/test/test_grudge.py b/test/test_grudge.py
index 6a42167ae71709d0a83980b7180645a8fcc07c22..ce0b199b8eae10688b77e91b2e08f50b4d0e48ea 100644
--- a/test/test_grudge.py
+++ b/test/test_grudge.py
@@ -31,8 +31,6 @@ from arraycontext import pytest_generate_tests_for_array_contexts
 pytest_generate_tests = pytest_generate_tests_for_array_contexts(
         [PytestPyOpenCLArrayContextFactory])
 
-from arraycontext.container.traversal import thaw
-
 from meshmode import _acf           # noqa: F401
 
 from meshmode.dof_array import flat_norm
@@ -92,12 +90,12 @@ def test_mass_mat_trig(actx_factory, ambient_dim, discr_tag):
         return actx.np.sin(x[0])**2
 
     volm_disc = dcoll.discr_from_dd(dof_desc.DD_VOLUME)
-    x_volm = thaw(volm_disc.nodes(), actx)
+    x_volm = actx.thaw(volm_disc.nodes())
     f_volm = f(x_volm)
     ones_volm = volm_disc.zeros(actx) + 1
 
     quad_disc = dcoll.discr_from_dd(dd_quad)
-    x_quad = thaw(quad_disc.nodes(), actx)
+    x_quad = actx.thaw(quad_disc.nodes())
     f_quad = f(x_quad)
     ones_quad = quad_disc.zeros(actx) + 1
 
@@ -271,7 +269,7 @@ def test_mass_operator_inverse(actx_factory, name):
             return actx.np.cos(4.0 * x[0])
 
         dd = dof_desc.DD_VOLUME
-        x_volm = thaw(volume_discr.nodes(), actx)
+        x_volm = actx.thaw(volume_discr.nodes())
         f_volm = f(x_volm)
         f_inv = op.inverse_mass(
             dcoll, op.mass(dcoll, dd, f_volm)
@@ -342,7 +340,7 @@ def test_face_normal_surface(actx_factory, mesh_name):
     )
     surf_normal = surf_normal / actx.np.sqrt(sum(surf_normal**2))
 
-    face_normal_i = thaw(dcoll.normal(df), actx)
+    face_normal_i = actx.thaw(dcoll.normal(df))
     face_normal_e = dcoll.opposite_face_connection()(face_normal_i)
 
     if mesh.ambient_dim == 3:
@@ -408,7 +406,7 @@ def test_tri_diff_mat(actx_factory, dim, order=4):
 
         dcoll = DiscretizationCollection(actx, mesh, order=4)
         volume_discr = dcoll.discr_from_dd(dof_desc.DD_VOLUME)
-        x = thaw(volume_discr.nodes(), actx)
+        x = actx.thaw(volume_discr.nodes())
 
         for axis in range(dim):
             df_num = op.local_grad(dcoll, f(x, axis))[axis]
@@ -450,7 +448,7 @@ def test_2d_gauss_theorem(actx_factory):
 
     dcoll = DiscretizationCollection(actx, mesh, order=2)
     volm_disc = dcoll.discr_from_dd(dof_desc.DD_VOLUME)
-    x_volm = thaw(volm_disc.nodes(), actx)
+    x_volm = actx.thaw(volm_disc.nodes())
 
     def f(x):
         return flat_obj_array(
@@ -462,7 +460,7 @@ def test_2d_gauss_theorem(actx_factory):
     int_1 = op.integral(dcoll, "vol", op.local_div(dcoll, f_volm))
 
     prj_f = op.project(dcoll, "vol", BTAG_ALL, f_volm)
-    normal = thaw(dcoll.normal(BTAG_ALL), actx)
+    normal = actx.thaw(dcoll.normal(BTAG_ALL))
     int_2 = op.integral(dcoll, BTAG_ALL, prj_f.dot(normal))
 
     assert abs(int_1 - int_2) < 1e-13
@@ -564,14 +562,14 @@ def test_surface_divergence_theorem(actx_factory, mesh_name, visualize=False):
         ambient_dim = dcoll.ambient_dim
 
         # variables
-        f_num = f(thaw(dcoll.nodes(dd=dd), actx))
-        f_quad_num = f(thaw(dcoll.nodes(dd=dq), actx))
+        f_num = f(actx.thaw(dcoll.nodes(dd=dd)))
+        f_quad_num = f(actx.thaw(dcoll.nodes(dd=dq)))
 
         from grudge.geometry import normal, summed_curvature
 
         kappa = summed_curvature(actx, dcoll, dd=dq)
         normal = normal(actx, dcoll, dd=dq)
-        face_normal = thaw(dcoll.normal(df), actx)
+        face_normal = actx.thaw(dcoll.normal(df))
         face_f = op.project(dcoll, dd, df, f_num)
 
         # operators
@@ -713,12 +711,12 @@ def test_convergence_advec(actx_factory, mesh_name, mesh_pars, op_type, flux_typ
                     "weak": WeakAdvectionOperator}[op_type]
         adv_operator = op_class(dcoll, v,
                                 inflow_u=lambda t: u_analytic(
-                                    thaw(dcoll.nodes(dd=BTAG_ALL), actx),
+                                    actx.thaw(dcoll.nodes(dd=BTAG_ALL)),
                                     t=t
                                 ),
                                 flux_type=flux_type)
 
-        nodes = thaw(dcoll.nodes(), actx)
+        nodes = actx.thaw(dcoll.nodes())
         u = u_analytic(nodes, t=0)
 
         def rhs(t, u):
@@ -812,7 +810,7 @@ def test_convergence_maxwell(actx_factory,  order):
         def analytic_sol(x, t=0):
             return get_rectangular_cavity_mode(actx, x, t, 1, (1, 2, 2))
 
-        nodes = thaw(dcoll.nodes(), actx)
+        nodes = actx.thaw(dcoll.nodes())
         fields = analytic_sol(nodes, t=0)
 
         from grudge.models.em import MaxwellOperator
@@ -909,7 +907,7 @@ def test_improvement_quadrature(actx_factory, order):
                 discr_tag_to_group_factory=discr_tag_to_group_factory
             )
 
-            nodes = thaw(dcoll.nodes(), actx)
+            nodes = actx.thaw(dcoll.nodes())
 
             def zero_inflow(dtag, t=0):
                 dd = dof_desc.DOFDesc(dtag, qtag)
@@ -959,7 +957,7 @@ def test_bessel(actx_factory):
 
     dcoll = DiscretizationCollection(actx, mesh, order=3)
 
-    nodes = thaw(dcoll.nodes(), actx)
+    nodes = actx.thaw(dcoll.nodes())
     r = actx.np.sqrt(nodes[0]**2 + nodes[1]**2)
 
     # FIXME: Bessel functions need to brought out of the symbolic
@@ -990,7 +988,7 @@ def test_norm_real(actx_factory, p):
             a=(0,)*dim, b=(1,)*dim,
             nelements_per_axis=(8,)*dim, order=1)
     dcoll = DiscretizationCollection(actx, mesh, order=4)
-    nodes = thaw(dcoll.nodes(), actx)
+    nodes = actx.thaw(dcoll.nodes())
 
     norm = op.norm(dcoll, nodes[0], p)
     if p == 2:
@@ -1011,7 +1009,7 @@ def test_norm_complex(actx_factory, p):
             a=(0,)*dim, b=(1,)*dim,
             nelements_per_axis=(8,)*dim, order=1)
     dcoll = DiscretizationCollection(actx, mesh, order=4)
-    nodes = thaw(dcoll.nodes(), actx)
+    nodes = actx.thaw(dcoll.nodes())
 
     norm = op.norm(dcoll, (1 + 1j)*nodes[0], p)
     if p == 2:
@@ -1032,7 +1030,7 @@ def test_norm_obj_array(actx_factory, p):
             a=(0,)*dim, b=(1,)*dim,
             nelements_per_axis=(8,)*dim, order=1)
     dcoll = DiscretizationCollection(actx, mesh, order=4)
-    nodes = thaw(dcoll.nodes(), actx)
+    nodes = actx.thaw(dcoll.nodes())
 
     norm = op.norm(dcoll, nodes, p)
 
diff --git a/test/test_modal_connections.py b/test/test_modal_connections.py
index 4f4354b2604d9bfb41c43cb233e5be5c694327f3..b633036f34050a2b25176d1676c000d85adb5096 100644
--- a/test/test_modal_connections.py
+++ b/test/test_modal_connections.py
@@ -25,7 +25,6 @@ from grudge.array_context import PytestPyOpenCLArrayContextFactory
 from arraycontext import pytest_generate_tests_for_array_contexts
 pytest_generate_tests = pytest_generate_tests_for_array_contexts(
         [PytestPyOpenCLArrayContextFactory])
-from arraycontext import thaw
 
 from meshmode.discretization.poly_element import (
     # Simplex group factories
@@ -76,7 +75,7 @@ def test_inverse_modal_connections(actx_factory, nodal_group_factory):
     dd_modal = dof_desc.DD_VOLUME_MODAL
     dd_volume = dof_desc.DD_VOLUME
 
-    x_nodal = thaw(dcoll.discr_from_dd(dd_volume).nodes()[0], actx)
+    x_nodal = actx.thaw(dcoll.discr_from_dd(dd_volume).nodes()[0])
     nodal_f = f(x_nodal)
 
     # Map nodal coefficients of f to modal coefficients
@@ -120,7 +119,7 @@ def test_inverse_modal_connections_quadgrid(actx_factory):
     dd_quad = dof_desc.DOFDesc(dof_desc.DTAG_VOLUME_ALL,
                                dof_desc.DISCR_TAG_QUAD)
 
-    x_quad = thaw(dcoll.discr_from_dd(dd_quad).nodes()[0], actx)
+    x_quad = actx.thaw(dcoll.discr_from_dd(dd_quad).nodes()[0])
     quad_f = f(x_quad)
 
     # Map nodal coefficients of f to modal coefficients
diff --git a/test/test_mpi_communication.py b/test/test_mpi_communication.py
index 8e053f6fda10314d569289fed7022bd428e4444d..47100b415f3f113de1d1b11e9d5b932af1fc5a8f 100644
--- a/test/test_mpi_communication.py
+++ b/test/test_mpi_communication.py
@@ -32,7 +32,6 @@ import logging
 import sys
 
 from grudge.array_context import MPIPyOpenCLArrayContext, MPIPytatoArrayContext
-from arraycontext.container.traversal import thaw, freeze
 
 logger = logging.getLogger(__name__)
 logging.basicConfig()
@@ -136,7 +135,7 @@ def _test_func_comparison_mpi_communication_entrypoint(actx):
 
     dcoll = DiscretizationCollection(actx, local_mesh, order=5)
 
-    x = thaw(dcoll.nodes(), actx)
+    x = actx.thaw(dcoll.nodes())
     myfunc = actx.np.sin(np.dot(x, [2, 3]))
 
     from grudge.dof_desc import as_dofdesc
@@ -214,7 +213,7 @@ def _test_mpi_wave_op_entrypoint(actx, visualize=False):
         source_center = np.array([0.1, 0.22, 0.33])[:dcoll.dim]
         source_width = 0.05
         source_omega = 3
-        nodes = thaw(dcoll.nodes(), actx)
+        nodes = actx.thaw(dcoll.nodes())
         source_center_dist = flat_obj_array(
             [nodes[i] - source_center[i] for i in range(dcoll.dim)]
         )
@@ -282,7 +281,7 @@ def _test_mpi_wave_op_entrypoint(actx, visualize=False):
     for step in range(nsteps):
         t = step*dt
         fields = rk4_step(fields, t=t, h=dt, f=compiled_rhs)
-        fields = thaw(freeze(fields, actx), actx)
+        fields = actx.thaw(actx.freeze(fields))
 
         norm = actx.to_numpy(op.norm(dcoll, fields, 2))
         logger.info("[%04d] t = %.5e |u| = %.5e elapsed %.5e",
diff --git a/test/test_op.py b/test/test_op.py
index 364a988f5e5107b188616917f5c233028267abb0..0eaac4f0e036168da2c92c0bb59a00a106089aa4 100644
--- a/test/test_op.py
+++ b/test/test_op.py
@@ -37,8 +37,6 @@ from arraycontext import pytest_generate_tests_for_array_contexts
 pytest_generate_tests = pytest_generate_tests_for_array_contexts(
         [PytestPyOpenCLArrayContextFactory])
 
-from arraycontext.container.traversal import thaw
-
 import logging
 
 logger = logging.getLogger(__name__)
@@ -89,7 +87,7 @@ def test_gradient(actx_factory, form, dim, order, vectorize, nested,
             result[dim-1] = result[dim-1] * (-np.pi/2*actx.np.sin(np.pi/2*x[dim-1]))
             return result
 
-        x = thaw(dcoll.nodes(), actx)
+        x = actx.thaw(dcoll.nodes())
 
         if vectorize:
             u = make_obj_array([(i+1)*f(x) for i in range(dim)])
@@ -99,7 +97,7 @@ def test_gradient(actx_factory, form, dim, order, vectorize, nested,
         def get_flux(u_tpair):
             dd = u_tpair.dd
             dd_allfaces = dd.with_dtag("all_faces")
-            normal = thaw(dcoll.normal(dd), actx)
+            normal = actx.thaw(dcoll.normal(dd))
             u_avg = u_tpair.avg
             if vectorize:
                 if nested:
@@ -213,7 +211,7 @@ def test_divergence(actx_factory, form, dim, order, vectorize, nested,
             result = result + deriv
             return result
 
-        x = thaw(dcoll.nodes(), actx)
+        x = actx.thaw(dcoll.nodes())
 
         if vectorize:
             u = make_obj_array([(i+1)*f(x) for i in range(dim)])
@@ -225,7 +223,7 @@ def test_divergence(actx_factory, form, dim, order, vectorize, nested,
         def get_flux(u_tpair):
             dd = u_tpair.dd
             dd_allfaces = dd.with_dtag("all_faces")
-            normal = thaw(dcoll.normal(dd), actx)
+            normal = actx.thaw(dcoll.normal(dd))
             flux = u_tpair.avg @ normal
             return op.project(dcoll, dd, dd_allfaces, flux)
 
diff --git a/test/test_reductions.py b/test/test_reductions.py
index d7f2f0b8c1ac413ae13eb7b73157d437735a83b8..9e7387badb5bafc31b0c4b10f7687e34919dc318 100644
--- a/test/test_reductions.py
+++ b/test/test_reductions.py
@@ -27,7 +27,6 @@ import numpy as np
 from dataclasses import dataclass
 
 from arraycontext import (
-    thaw,
     with_container_arithmetic,
     dataclass_array_container,
     pytest_generate_tests_for_array_contexts
@@ -70,7 +69,7 @@ def test_nodal_reductions(actx_factory, mesh_size, with_initial):
 
     mesh = builder.get_mesh(mesh_size, builder.mesh_order)
     dcoll = DiscretizationCollection(actx, mesh, order=builder.order)
-    x = thaw(dcoll.nodes(), actx)
+    x = actx.thaw(dcoll.nodes())
 
     def f(x):
         return -actx.np.sin(10*x[0])
@@ -137,7 +136,7 @@ def test_elementwise_reductions(actx_factory):
     nelements = 4
     mesh = builder.get_mesh(nelements, builder.mesh_order)
     dcoll = DiscretizationCollection(actx, mesh, order=builder.order)
-    x = thaw(dcoll.nodes(), actx)
+    x = actx.thaw(dcoll.nodes())
 
     def f(x):
         return actx.np.sin(x[0])
@@ -197,7 +196,7 @@ def test_nodal_reductions_with_container(actx_factory):
 
     mesh = builder.get_mesh(4, builder.mesh_order)
     dcoll = DiscretizationCollection(actx, mesh, order=builder.order)
-    x = thaw(dcoll.nodes(), actx)
+    x = actx.thaw(dcoll.nodes())
 
     def f(x):
         return -actx.np.sin(10*x[0]) * actx.np.cos(2*x[1])
@@ -245,7 +244,7 @@ def test_elementwise_reductions_with_container(actx_factory):
     nelements = 4
     mesh = builder.get_mesh(nelements, builder.mesh_order)
     dcoll = DiscretizationCollection(actx, mesh, order=builder.order)
-    x = thaw(dcoll.nodes(), actx)
+    x = actx.thaw(dcoll.nodes())
 
     def f(x):
         return actx.np.sin(x[0]) * actx.np.sin(x[1])