diff --git a/examples/advection/var-velocity.py b/examples/advection/var-velocity.py
index 454453972795760b674b25cbcd0fdf417bb01ed5..388c4ce576fd8e248432831f14efaadb248fb4a1 100644
--- a/examples/advection/var-velocity.py
+++ b/examples/advection/var-velocity.py
@@ -26,8 +26,8 @@ import os
 import numpy as np
 
 import pyopencl as cl
-import pyopencl.array
-import pyopencl.clmath
+
+from meshmode.array_context import PyOpenCLArrayContext
 
 from grudge import bind, sym
 from pytools.obj_array import flat_obj_array
@@ -92,6 +92,7 @@ class Plotter:
 def main(ctx_factory, dim=2, order=4, product_tag=None, visualize=True):
     cl_ctx = ctx_factory()
     queue = cl.CommandQueue(cl_ctx)
+    actx = PyOpenCLArrayContext(queue)
 
     # {{{ parameters
 
@@ -145,7 +146,7 @@ def main(ctx_factory, dim=2, order=4, product_tag=None, visualize=True):
         quad_tag_to_group_factory = {}
 
     from grudge import DGDiscretizationWithBoundaries
-    discr = DGDiscretizationWithBoundaries(cl_ctx, mesh, order=order,
+    discr = DGDiscretizationWithBoundaries(actx, mesh, order=order,
             quad_tag_to_group_factory=quad_tag_to_group_factory)
 
     # }}}
@@ -176,10 +177,10 @@ def main(ctx_factory, dim=2, order=4, product_tag=None, visualize=True):
             flux_type=flux_type)
 
     bound_op = bind(discr, op.sym_operator())
-    u = bind(discr, f_gaussian(sym.nodes(dim)))(queue, t=0)
+    u = bind(discr, f_gaussian(sym.nodes(dim)))(actx, t=0)
 
     def rhs(t, u):
-        return bound_op(queue, t=t, u=u)
+        return bound_op(t=t, u=u)
 
     # }}}
 
@@ -197,7 +198,7 @@ def main(ctx_factory, dim=2, order=4, product_tag=None, visualize=True):
             continue
 
         if step % 10 == 0:
-            norm_u = norm(queue, u=event.state_component)
+            norm_u = norm(u=event.state_component)
             plot(event, "fld-var-velocity-%04d" % step)
 
         step += 1
diff --git a/examples/advection/weak.py b/examples/advection/weak.py
index 14ed08e1a757d149c1c96ecb9d84eb60d15400b1..88703a0e930ff703d15e3e9dee6d7be11974002a 100644
--- a/examples/advection/weak.py
+++ b/examples/advection/weak.py
@@ -22,8 +22,8 @@ import numpy as np
 import numpy.linalg as la
 
 import pyopencl as cl
-import pyopencl.array
-import pyopencl.clmath
+
+from meshmode.array_context import PyOpenCLArrayContext
 
 from grudge import bind, sym
 
@@ -88,6 +88,7 @@ class Plotter:
 def main(ctx_factory, dim=2, order=4, visualize=True):
     cl_ctx = ctx_factory()
     queue = cl.CommandQueue(cl_ctx)
+    actx = PyOpenCLArrayContext(queue)
 
     # {{{ parameters
 
@@ -123,7 +124,7 @@ def main(ctx_factory, dim=2, order=4, visualize=True):
             order=order)
 
     from grudge import DGDiscretizationWithBoundaries
-    discr = DGDiscretizationWithBoundaries(cl_ctx, mesh, order=order)
+    discr = DGDiscretizationWithBoundaries(actx, mesh, order=order)
 
     # }}}
 
@@ -142,10 +143,10 @@ def main(ctx_factory, dim=2, order=4, visualize=True):
         flux_type=flux_type)
 
     bound_op = bind(discr, op.sym_operator())
-    u = bind(discr, u_analytic(sym.nodes(dim)))(queue, t=0)
+    u = bind(discr, u_analytic(sym.nodes(dim)))(actx, t=0)
 
     def rhs(t, u):
-        return bound_op(queue, t=t, u=u)
+        return bound_op(t=t, u=u)
 
     # }}}
 
@@ -165,7 +166,7 @@ def main(ctx_factory, dim=2, order=4, visualize=True):
             continue
 
         if step % 10 == 0:
-            norm_u = norm(queue, u=event.state_component)
+            norm_u = norm(u=event.state_component)
             plot(event, "fld-weak-%04d" % step)
 
         step += 1
diff --git a/examples/geometry.py b/examples/geometry.py
index 6e146f21dd31545b392d8006d45d94bed7a9db02..df81298d4865b0dce384952633ee4fd00717c72b 100644
--- a/examples/geometry.py
+++ b/examples/geometry.py
@@ -28,15 +28,18 @@ import numpy as np  # noqa
 import pyopencl as cl
 from grudge import sym, bind, DGDiscretizationWithBoundaries, shortcuts
 
+from meshmode.array_context import PyOpenCLArrayContext
+
 
 def main(write_output=True):
     cl_ctx = cl.create_some_context()
     queue = cl.CommandQueue(cl_ctx)
+    actx = PyOpenCLArrayContext(queue)
 
     from meshmode.mesh.generation import generate_warped_rect_mesh
     mesh = generate_warped_rect_mesh(dim=2, order=4, n=6)
 
-    discr = DGDiscretizationWithBoundaries(cl_ctx, mesh, order=4)
+    discr = DGDiscretizationWithBoundaries(actx, mesh, order=4)
 
     sym_op = sym.normal(sym.BTAG_ALL, mesh.dim)
     #sym_op = sym.nodes(mesh.dim, where=sym.BTAG_ALL)
@@ -45,7 +48,7 @@ def main(write_output=True):
     print()
     print(op.eval_code)
 
-    vec = op(queue)
+    vec = op(actx)
 
     vis = shortcuts.make_visualizer(discr, 4)
     vis.write_vtk_file("geo.vtu", [
diff --git a/examples/maxwell/cavities.py b/examples/maxwell/cavities.py
index a58f27399ccfb6e32a16d6cf594ed2a3c78076b2..eb2ee28dae6b2ed74c2760d2359d9542a9df2f59 100644
--- a/examples/maxwell/cavities.py
+++ b/examples/maxwell/cavities.py
@@ -27,6 +27,9 @@ THE SOFTWARE.
 
 import numpy as np
 import pyopencl as cl
+
+from meshmode.array_context import PyOpenCLArrayContext
+
 from grudge.shortcuts import set_up_rk4
 from grudge import sym, bind, DGDiscretizationWithBoundaries
 
@@ -39,6 +42,7 @@ STEPS = 60
 def main(dims, write_output=True, order=4):
     cl_ctx = cl.create_some_context()
     queue = cl.CommandQueue(cl_ctx)
+    actx = PyOpenCLArrayContext(queue)
 
     from meshmode.mesh.generation import generate_regular_rect_mesh
     mesh = generate_regular_rect_mesh(
@@ -46,7 +50,7 @@ def main(dims, write_output=True, order=4):
             b=(1.0,)*dims,
             n=(5,)*dims)
 
-    discr = DGDiscretizationWithBoundaries(cl_ctx, mesh, order=order)
+    discr = DGDiscretizationWithBoundaries(actx, mesh, order=order)
 
     if 0:
         epsilon0 = 8.8541878176e-12  # C**2 / (N m**2)
@@ -60,14 +64,12 @@ def main(dims, write_output=True, order=4):
     from grudge.models.em import MaxwellOperator
     op = MaxwellOperator(epsilon, mu, flux_type=0.5, dimensions=dims)
 
-    queue = cl.CommandQueue(discr.cl_context)
-
     if dims == 3:
         sym_mode = get_rectangular_cavity_mode(1, (1, 2, 2))
-        fields = bind(discr, sym_mode)(queue, t=0, epsilon=epsilon, mu=mu)
+        fields = bind(discr, sym_mode)(actx, t=0, epsilon=epsilon, mu=mu)
     else:
         sym_mode = get_rectangular_cavity_mode(1, (2, 3))
-        fields = bind(discr, sym_mode)(queue, t=0)
+        fields = bind(discr, sym_mode)(actx, t=0)
 
     # FIXME
     #dt = op.estimate_rk4_timestep(discr, fields=fields)
@@ -78,7 +80,7 @@ def main(dims, write_output=True, order=4):
     bound_op = bind(discr, op.sym_operator())
 
     def rhs(t, w):
-        return bound_op(queue, t=t, w=w)
+        return bound_op(t=t, w=w)
 
     if mesh.dim == 2:
         dt = 0.004
@@ -117,8 +119,8 @@ def main(dims, write_output=True, order=4):
 
             step += 1
 
-            print(step, event.t, norm(queue, u=e[0]), norm(queue, u=e[1]),
-                    norm(queue, u=h[0]), norm(queue, u=h[1]),
+            print(step, event.t, norm(u=e[0]), norm(u=e[1]),
+                    norm(u=h[0]), norm(u=h[1]),
                     time()-t_last_step)
             if step % 10 == 0:
                 e, h = op.split_eh(event.state_component)
diff --git a/examples/wave/var-propagation-speed.py b/examples/wave/var-propagation-speed.py
index 153b85d2f1b4d822a0d16b3e6c6eec74f8bdf1ce..f392b420d4e06c6fb547c6461c0290726388cf29 100644
--- a/examples/wave/var-propagation-speed.py
+++ b/examples/wave/var-propagation-speed.py
@@ -30,10 +30,13 @@ import pyopencl as cl
 from grudge.shortcuts import set_up_rk4
 from grudge import sym, bind, DGDiscretizationWithBoundaries
 
+from meshmode.array_context import PyOpenCLArrayContext
+
 
 def main(write_output=True, order=4):
     cl_ctx = cl.create_some_context()
     queue = cl.CommandQueue(cl_ctx)
+    actx = PyOpenCLArrayContext(queue)
 
     dims = 2
     from meshmode.mesh.generation import generate_regular_rect_mesh
@@ -42,7 +45,7 @@ def main(write_output=True, order=4):
             b=(0.5,)*dims,
             n=(20,)*dims)
 
-    discr = DGDiscretizationWithBoundaries(cl_ctx, mesh, order=order)
+    discr = DGDiscretizationWithBoundaries(actx, mesh, order=order)
 
     source_center = np.array([0.1, 0.22, 0.33])[:mesh.dim]
     source_width = 0.05
@@ -69,19 +72,18 @@ def main(write_output=True, order=4):
             radiation_tag=BTAG_ALL,
             flux_type="upwind")
 
-    queue = cl.CommandQueue(discr.cl_context)
     from pytools.obj_array import flat_obj_array
-    fields = flat_obj_array(discr.zeros(queue),
-            [discr.zeros(queue) for i in range(discr.dim)])
+    fields = flat_obj_array(discr.zeros(actx),
+            [discr.zeros(actx) for i in range(discr.dim)])
 
     op.check_bc_coverage(mesh)
 
-    c_eval = bind(discr, c)(queue)
+    c_eval = bind(discr, c)(actx)
 
     bound_op = bind(discr, op.sym_operator())
 
     def rhs(t, w):
-        return bound_op(queue, t=t, w=w)
+        return bound_op(t=t, w=w)
 
     if mesh.dim == 2:
         dt = 0.04 * 0.3
@@ -110,7 +112,7 @@ def main(write_output=True, order=4):
 
             step += 1
 
-            print(step, event.t, norm(queue, u=event.state_component[0]),
+            print(step, event.t, norm(u=event.state_component[0]),
                     time()-t_last_step)
             if step % 10 == 0:
                 vis.write_vtk_file("fld-var-propogation-speed-%04d.vtu" % step,
diff --git a/examples/wave/wave-eager.py b/examples/wave/wave-eager.py
index 5737dcd57ed41482ecaf96749e7923ae235d0b62..76dab5095d2f0f53e8bf9bace9a02b4498165f62 100644
--- a/examples/wave/wave-eager.py
+++ b/examples/wave/wave-eager.py
@@ -26,23 +26,17 @@ THE SOFTWARE.
 import numpy as np
 import numpy.linalg as la  # noqa
 import pyopencl as cl
-import pyopencl.array as cla  # noqa
-import pyopencl.clmath as clmath
-from pytools.obj_array import (
-        with_object_array_or_scalar)
-from meshmode.mesh import BTAG_ALL, BTAG_NONE  # noqa
-from grudge.eager import EagerDGDiscretization, with_queue
-from grudge.symbolic.primitives import TracePair
-from grudge.shortcuts import make_visualizer
 
 from pytools.obj_array import flat_obj_array, make_obj_array
 
-def interior_trace_pair(discr, vec):
-    i = discr.interp("vol", "int_faces", vec)
-    e = with_object_array_or_scalar(
-            lambda el: discr.opposite_face_connection()(el.queue, el),
-            i)
-    return TracePair("int_faces", i, e)
+from meshmode.array_context import PyOpenCLArrayContext
+from meshmode.dof_array import thaw
+
+from meshmode.mesh import BTAG_ALL, BTAG_NONE  # noqa
+
+from grudge.eager import EagerDGDiscretization, interior_trace_pair
+from grudge.shortcuts import make_visualizer
+from grudge.symbolic.primitives import TracePair
 
 
 # {{{ wave equation bits
@@ -51,7 +45,7 @@ def wave_flux(discr, c, w_tpair):
     u = w_tpair[0]
     v = w_tpair[1:]
 
-    normal = with_queue(u.int.queue, discr.normal(w_tpair.dd))
+    normal = thaw(u.int.array_context, discr.normal(w_tpair.dd))
 
     def normal_times(scalar):
         # workaround for object array behavior
@@ -106,20 +100,21 @@ def rk4_step(y, t, h, f):
     return y + h/6*(k1 + 2*k2 + 2*k3 + k4)
 
 
-def bump(discr, queue, t=0):
+def bump(discr, actx, t=0):
     source_center = np.array([0.2, 0.35, 0.1])[:discr.dim]
     source_width = 0.05
     source_omega = 3
 
-    nodes = discr.nodes().with_queue(queue)
+    nodes = thaw(actx, discr.nodes())
     center_dist = flat_obj_array([
         nodes[i] - source_center[i]
         for i in range(discr.dim)
         ])
 
+    exp = actx.special_func("exp")
     return (
         np.cos(source_omega*t)
-        * clmath.exp(
+        * exp(
             -np.dot(center_dist, center_dist)
             / source_width**2))
 
@@ -127,6 +122,7 @@ def bump(discr, queue, t=0):
 def main():
     cl_ctx = cl.create_some_context()
     queue = cl.CommandQueue(cl_ctx)
+    actx = PyOpenCLArrayContext(queue)
 
     dim = 2
     nel_1d = 16
@@ -149,11 +145,11 @@ def main():
 
     print("%d elements" % mesh.nelements)
 
-    discr = EagerDGDiscretization(cl_ctx, mesh, order=order)
+    discr = EagerDGDiscretization(actx, mesh, order=order)
 
     fields = flat_obj_array(
-            bump(discr, queue),
-            [discr.zeros(queue) for i in range(discr.dim)]
+            bump(discr, actx),
+            [discr.zeros(actx) for i in range(discr.dim)]
             )
 
     vis = make_visualizer(discr, discr.order+3 if dim == 2 else discr.order)
@@ -168,7 +164,7 @@ def main():
         fields = rk4_step(fields, t, dt, rhs)
 
         if istep % 10 == 0:
-            print(istep, t, la.norm(fields[0].get()))
+            print(istep, t, discr.norm(fields[0]))
             vis.write_vtk_file("fld-wave-eager-%04d.vtu" % istep,
                     [
                         ("u", fields[0]),
diff --git a/examples/wave/wave-min-mpi.py b/examples/wave/wave-min-mpi.py
index 9a8e70bf8b2f2ee76bca2c586089f7de30958f5a..bca2610199ed1cf50ae0c178d28caf1579730e3c 100644
--- a/examples/wave/wave-min-mpi.py
+++ b/examples/wave/wave-min-mpi.py
@@ -27,6 +27,7 @@ THE SOFTWARE.
 
 import numpy as np
 import pyopencl as cl
+from meshmode.array_context import PyOpenCLArrayContext
 from grudge.shortcuts import set_up_rk4
 from grudge import sym, bind, DGDiscretizationWithBoundaries
 from mpi4py import MPI
@@ -35,6 +36,7 @@ from mpi4py import MPI
 def main(write_output=True, order=4):
     cl_ctx = cl.create_some_context()
     queue = cl.CommandQueue(cl_ctx)
+    actx = PyOpenCLArrayContext(queue)
 
     comm = MPI.COMM_WORLD
     num_parts = comm.Get_size()
@@ -61,7 +63,7 @@ def main(write_output=True, order=4):
     else:
         local_mesh = mesh_dist.receive_mesh_part()
 
-    discr = DGDiscretizationWithBoundaries(cl_ctx, local_mesh, order=order,
+    discr = DGDiscretizationWithBoundaries(actx, local_mesh, order=order,
             mpi_communicator=comm)
 
     if local_mesh.dim == 2:
@@ -90,10 +92,10 @@ def main(write_output=True, order=4):
             radiation_tag=BTAG_ALL,
             flux_type="upwind")
 
-    queue = cl.CommandQueue(discr.cl_context)
     from pytools.obj_array import flat_obj_array
-    fields = flat_obj_array(discr.zeros(queue),
-            [discr.zeros(queue) for i in range(discr.dim)])
+    fields = flat_obj_array(
+            discr.zeros(actx),
+            [discr.zeros(actx) for i in range(discr.dim)])
 
     # FIXME
     #dt = op.estimate_rk4_timestep(discr, fields=fields)
@@ -104,7 +106,7 @@ def main(write_output=True, order=4):
     bound_op = bind(discr, op.sym_operator())
 
     def rhs(t, w):
-        return bound_op(queue, t=t, w=w)
+        return bound_op(t=t, w=w)
 
     dt_stepper = set_up_rk4("w", dt, fields, rhs)
 
@@ -130,7 +132,7 @@ def main(write_output=True, order=4):
 
             step += 1
 
-            print(step, event.t, norm(queue, u=event.state_component[0]),
+            print(step, event.t, norm(u=event.state_component[0]),
                     time()-t_last_step)
             if step % 10 == 0:
                 vis.write_vtk_file(
diff --git a/examples/wave/wave-min.py b/examples/wave/wave-min.py
index 665ce71f254a61f105b6220d0eff07d91a6b2429..de6c9f92e75d449e6a42de80f7b6fa7270df8a15 100644
--- a/examples/wave/wave-min.py
+++ b/examples/wave/wave-min.py
@@ -27,6 +27,7 @@ THE SOFTWARE.
 
 import numpy as np
 import pyopencl as cl
+from meshmode.array_context import PyOpenCLArrayContext
 from grudge.shortcuts import set_up_rk4
 from grudge import sym, bind, DGDiscretizationWithBoundaries
 
@@ -34,6 +35,7 @@ from grudge import sym, bind, DGDiscretizationWithBoundaries
 def main(write_output=True, order=4):
     cl_ctx = cl.create_some_context()
     queue = cl.CommandQueue(cl_ctx)
+    actx = PyOpenCLArrayContext(queue)
 
     dims = 2
     from meshmode.mesh.generation import generate_regular_rect_mesh
@@ -49,7 +51,7 @@ def main(write_output=True, order=4):
 
     print("%d elements" % mesh.nelements)
 
-    discr = DGDiscretizationWithBoundaries(cl_ctx, mesh, order=order)
+    discr = DGDiscretizationWithBoundaries(actx, mesh, order=order)
 
     source_center = np.array([0.1, 0.22, 0.33])[:mesh.dim]
     source_width = 0.05
@@ -72,10 +74,9 @@ def main(write_output=True, order=4):
             radiation_tag=BTAG_ALL,
             flux_type="upwind")
 
-    queue = cl.CommandQueue(discr.cl_context)
     from pytools.obj_array import flat_obj_array
-    fields = flat_obj_array(discr.zeros(queue),
-            [discr.zeros(queue) for i in range(discr.dim)])
+    fields = flat_obj_array(discr.zeros(actx),
+            [discr.zeros(actx) for i in range(discr.dim)])
 
     # FIXME
     #dt = op.estimate_rk4_timestep(discr, fields=fields)
@@ -86,7 +87,7 @@ def main(write_output=True, order=4):
     bound_op = bind(discr, op.sym_operator())
 
     def rhs(t, w):
-        return bound_op(queue, t=t, w=w)
+        return bound_op(t=t, w=w)
 
     dt_stepper = set_up_rk4("w", dt, fields, rhs)
 
@@ -110,7 +111,7 @@ def main(write_output=True, order=4):
 
             step += 1
 
-            print(step, event.t, norm(queue, u=event.state_component[0]),
+            print(step, event.t, norm(u=event.state_component[0]),
                     time()-t_last_step)
             if step % 10 == 0:
                 vis.write_vtk_file("fld-wave-min-%04d.vtu" % step,