diff --git a/grudge/discretization.py b/grudge/discretization.py
index 64c00e94f7b3668b8a8f1ae96ac59a1b8f8daf66..a8bee7e11a046aaeaf810f3f47c0e01ae6df995d 100644
--- a/grudge/discretization.py
+++ b/grudge/discretization.py
@@ -25,9 +25,9 @@ THE SOFTWARE.
 
 import six
 from pytools import memoize_method
-import pyopencl as cl
 from grudge import sym
-import numpy as np
+import numpy as np  # noqa: F401
+from meshmode.array_context import ArrayContext
 
 
 # FIXME Naming not ideal
@@ -40,7 +40,6 @@ class DGDiscretizationWithBoundaries(DiscretizationBase):
     .. automethod :: discr_from_dd
     .. automethod :: connection_from_dds
 
-    .. autoattribute :: cl_context
     .. autoattribute :: dim
     .. autoattribute :: ambient_dim
     .. autoattribute :: mesh
@@ -49,7 +48,7 @@ class DGDiscretizationWithBoundaries(DiscretizationBase):
     .. automethod :: zeros
     """
 
-    def __init__(self, cl_ctx, mesh, order, quad_tag_to_group_factory=None,
+    def __init__(self, array_context, mesh, order, quad_tag_to_group_factory=None,
             mpi_communicator=None):
         """
         :param quad_tag_to_group_factory: A mapping from quadrature tags (typically
@@ -60,6 +59,8 @@ class DGDiscretizationWithBoundaries(DiscretizationBase):
             be carried out with the standard volume discretization.
         """
 
+        self._setup_actx = array_context
+
         if quad_tag_to_group_factory is None:
             quad_tag_to_group_factory = {}
 
@@ -68,7 +69,7 @@ class DGDiscretizationWithBoundaries(DiscretizationBase):
 
         from meshmode.discretization import Discretization
 
-        self._volume_discr = Discretization(cl_ctx, mesh,
+        self._volume_discr = Discretization(array_context, mesh,
                 self.group_factory_for_quadrature_tag(sym.QTAG_NONE))
 
         # {{{ management of discretization-scoped common subexpressions
@@ -81,9 +82,9 @@ class DGDiscretizationWithBoundaries(DiscretizationBase):
 
         # }}}
 
-        with cl.CommandQueue(cl_ctx) as queue:
-            self._dist_boundary_connections = \
-                    self._set_up_distributed_communication(mpi_communicator, queue)
+        self._dist_boundary_connections = \
+                self._set_up_distributed_communication(
+                        mpi_communicator, array_context)
 
         self.mpi_communicator = mpi_communicator
 
@@ -97,7 +98,7 @@ class DGDiscretizationWithBoundaries(DiscretizationBase):
             return self.mpi_communicator.Get_rank() \
                     == self._get_management_rank_index()
 
-    def _set_up_distributed_communication(self, mpi_communicator, queue):
+    def _set_up_distributed_communication(self, mpi_communicator, array_context):
         from_dd = sym.DOFDesc("vol", sym.QTAG_NONE)
 
         from meshmode.distributed import get_connected_partitions
@@ -118,7 +119,8 @@ class DGDiscretizationWithBoundaries(DiscretizationBase):
                     from_dd,
                     sym.DOFDesc(sym.BTAG_PARTITION(i_remote_part), sym.QTAG_NONE))
             setup_helper = setup_helpers[i_remote_part] = MPIBoundaryCommSetupHelper(
-                    mpi_communicator, queue, conn, i_remote_part, grp_factory)
+                    mpi_communicator, array_context, conn,
+                    i_remote_part, grp_factory)
             setup_helper.post_sends()
 
         for i_remote_part, setup_helper in six.iteritems(setup_helpers):
@@ -152,7 +154,7 @@ class DGDiscretizationWithBoundaries(DiscretizationBase):
 
             from meshmode.discretization import Discretization
             return Discretization(
-                    self._volume_discr.cl_context,
+                    self._setup_actx,
                     no_quad_discr.mesh,
                     self.group_factory_for_quadrature_tag(qtag))
 
@@ -186,6 +188,7 @@ class DGDiscretizationWithBoundaries(DiscretizationBase):
                     make_face_to_all_faces_embedding
 
             return make_face_to_all_faces_embedding(
+                    self._setup_actx,
                     faces_conn, self.discr_from_dd(to_dd),
                     self.discr_from_dd(from_dd))
 
@@ -224,6 +227,7 @@ class DGDiscretizationWithBoundaries(DiscretizationBase):
                     make_same_mesh_connection
 
             return make_same_mesh_connection(
+                    self._setup_actx,
                     self.discr_from_dd(to_dd),
                     self.discr_from_dd(from_dd))
 
@@ -276,7 +280,7 @@ class DGDiscretizationWithBoundaries(DiscretizationBase):
     def _quad_volume_discr(self, quadrature_tag):
         from meshmode.discretization import Discretization
 
-        return Discretization(self._volume_discr.cl_context, self._volume_discr.mesh,
+        return Discretization(self._setup_actx, self._volume_discr.mesh,
                 self.group_factory_for_quadrature_tag(quadrature_tag))
 
     # {{{ boundary
@@ -285,9 +289,10 @@ class DGDiscretizationWithBoundaries(DiscretizationBase):
     def _boundary_connection(self, boundary_tag):
         from meshmode.discretization.connection import make_face_restriction
         return make_face_restriction(
-                        self._volume_discr,
-                        self.group_factory_for_quadrature_tag(sym.QTAG_NONE),
-                        boundary_tag=boundary_tag)
+                self._setup_actx,
+                self._volume_discr,
+                self.group_factory_for_quadrature_tag(sym.QTAG_NONE),
+                boundary_tag=boundary_tag)
 
     # }}}
 
@@ -298,21 +303,24 @@ class DGDiscretizationWithBoundaries(DiscretizationBase):
         from meshmode.discretization.connection import (
                 make_face_restriction, FACE_RESTR_INTERIOR)
         return make_face_restriction(
-                        self._volume_discr,
-                        self.group_factory_for_quadrature_tag(sym.QTAG_NONE),
-                        FACE_RESTR_INTERIOR,
+                self._setup_actx,
+                self._volume_discr,
+                self.group_factory_for_quadrature_tag(sym.QTAG_NONE),
+                FACE_RESTR_INTERIOR,
 
-                        # FIXME: This will need to change as soon as we support
-                        # pyramids or other elements with non-identical face
-                        # types.
-                        per_face_groups=False)
+                # FIXME: This will need to change as soon as we support
+                # pyramids or other elements with non-identical face
+                # types.
+                per_face_groups=False)
 
     @memoize_method
     def opposite_face_connection(self):
         from meshmode.discretization.connection import \
                 make_opposite_face_connection
 
-        return make_opposite_face_connection(self._interior_faces_connection())
+        return make_opposite_face_connection(
+                self._setup_actx,
+                self._interior_faces_connection())
 
     # }}}
 
@@ -323,21 +331,18 @@ class DGDiscretizationWithBoundaries(DiscretizationBase):
         from meshmode.discretization.connection import (
                 make_face_restriction, FACE_RESTR_ALL)
         return make_face_restriction(
-                        self._volume_discr,
-                        self.group_factory_for_quadrature_tag(sym.QTAG_NONE),
-                        FACE_RESTR_ALL,
+                self._setup_actx,
+                self._volume_discr,
+                self.group_factory_for_quadrature_tag(sym.QTAG_NONE),
+                FACE_RESTR_ALL,
 
-                        # FIXME: This will need to change as soon as we support
-                        # pyramids or other elements with non-identical face
-                        # types.
-                        per_face_groups=False)
+                # FIXME: This will need to change as soon as we support
+                # pyramids or other elements with non-identical face
+                # types.
+                per_face_groups=False)
 
     # }}}
 
-    @property
-    def cl_context(self):
-        return self._volume_discr.cl_context
-
     @property
     def dim(self):
         return self._volume_discr.dim
@@ -358,13 +363,11 @@ class DGDiscretizationWithBoundaries(DiscretizationBase):
     def mesh(self):
         return self._volume_discr.mesh
 
-    def empty(self, queue=None, dtype=None, extra_dims=None, allocator=None):
-        return self._volume_discr.empty(queue, dtype, extra_dims=extra_dims,
-                allocator=allocator)
+    def empty(self, array_context: ArrayContext, dtype=None):
+        return self._volume_discr.empty(array_context, dtype)
 
-    def zeros(self, queue, dtype=None, extra_dims=None, allocator=None):
-        return self._volume_discr.zeros(queue, dtype, extra_dims=extra_dims,
-                allocator=allocator)
+    def zeros(self, array_context: ArrayContext, dtype=None):
+        return self._volume_discr.zeros(array_context, dtype)
 
     def is_volume_where(self, where):
         from grudge import sym
diff --git a/grudge/eager.py b/grudge/eager.py
index e64e2fa8344a848f4d0a218cfd0286ec9f896e2a..8b05edd0cebeb4541c520ae3a31ea69b292c3ea1 100644
--- a/grudge/eager.py
+++ b/grudge/eager.py
@@ -26,31 +26,23 @@ THE SOFTWARE.
 import numpy as np  # noqa
 from grudge.discretization import DGDiscretizationWithBoundaries
 from pytools import memoize_method
-from pytools.obj_array import (
-        with_object_array_or_scalar,
-        is_obj_array)
-import pyopencl as cl
+from pytools.obj_array import obj_array_vectorize
 import pyopencl.array as cla  # noqa
 from grudge import sym, bind
-from meshmode.mesh import BTAG_ALL, BTAG_NONE  # noqa
-
-
-def with_queue(queue, ary):
-    return with_object_array_or_scalar(
-            lambda x: x.with_queue(queue), ary)
-
 
-def without_queue(ary):
-    return with_queue(None, ary)
+from meshmode.mesh import BTAG_ALL, BTAG_NONE  # noqa
+from meshmode.dof_array import freeze, DOFArray
 
 
 class EagerDGDiscretization(DGDiscretizationWithBoundaries):
     def interp(self, src, tgt, vec):
-        if is_obj_array(vec):
-            return with_object_array_or_scalar(
+        if (isinstance(vec, np.ndarray)
+                and vec.dtype.char == "O"
+                and not isinstance(vec, DOFArray)):
+            return obj_array_vectorize(
                     lambda el: self.interp(src, tgt, el), vec)
 
-        return self.connection_from_dds(src, tgt)(vec.queue, vec)
+        return self.connection_from_dds(src, tgt)(vec)
 
     def nodes(self):
         return self._volume_discr.nodes()
@@ -60,7 +52,7 @@ class EagerDGDiscretization(DGDiscretizationWithBoundaries):
         return bind(self, sym.nabla(self.dim) * sym.Variable("u"))
 
     def grad(self, vec):
-        return self._bound_grad()(vec.queue, u=vec)
+        return self._bound_grad()(u=vec)
 
     def div(self, vecs):
         return sum(
@@ -71,7 +63,7 @@ class EagerDGDiscretization(DGDiscretizationWithBoundaries):
         return bind(self, sym.stiffness_t(self.dim) * sym.Variable("u"))
 
     def weak_grad(self, vec):
-        return self._bound_weak_grad()(vec.queue, u=vec)
+        return self._bound_weak_grad()(u=vec)
 
     def weak_div(self, vecs):
         return sum(
@@ -79,22 +71,25 @@ class EagerDGDiscretization(DGDiscretizationWithBoundaries):
 
     @memoize_method
     def normal(self, dd):
-        with cl.CommandQueue(self.cl_context) as queue:
-            surface_discr = self.discr_from_dd(dd)
-            return without_queue(
-                    bind(self, sym.normal(
-                        dd, surface_discr.ambient_dim, surface_discr.dim))(queue))
+        surface_discr = self.discr_from_dd(dd)
+        actx = surface_discr._setup_actx
+        return freeze(
+                bind(self,
+                    sym.normal(dd, surface_discr.ambient_dim, surface_discr.dim))
+                (array_context=actx))
 
     @memoize_method
     def _bound_inverse_mass(self):
         return bind(self, sym.InverseMassOperator()(sym.Variable("u")))
 
     def inverse_mass(self, vec):
-        if is_obj_array(vec):
-            return with_object_array_or_scalar(
+        if (isinstance(vec, np.ndarray)
+                and vec.dtype.char == "O"
+                and not isinstance(vec, DOFArray)):
+            return obj_array_vectorize(
                     lambda el: self.inverse_mass(el), vec)
 
-        return self._bound_inverse_mass()(vec.queue, u=vec)
+        return self._bound_inverse_mass()(u=vec)
 
     @memoize_method
     def _bound_face_mass(self):
@@ -102,10 +97,34 @@ class EagerDGDiscretization(DGDiscretizationWithBoundaries):
         return bind(self, sym.FaceMassOperator()(u))
 
     def face_mass(self, vec):
-        if is_obj_array(vec):
-            return with_object_array_or_scalar(
+        if (isinstance(vec, np.ndarray)
+                and vec.dtype.char == "O"
+                and not isinstance(vec, DOFArray)):
+            return obj_array_vectorize(
                     lambda el: self.face_mass(el), vec)
 
-        return self._bound_face_mass()(vec.queue, u=vec)
+        return self._bound_face_mass()(u=vec)
+
+    @memoize_method
+    def _norm(self, p):
+        return bind(self, sym.norm(p, sym.var("arg")))
+
+    def norm(self, vec, p=2):
+        return self._norm(p)(arg=vec)
+
+
+def interior_trace_pair(discr, vec):
+    i = discr.interp("vol", "int_faces", vec)
+
+    if (isinstance(vec, np.ndarray)
+            and vec.dtype.char == "O"
+            and not isinstance(vec, DOFArray)):
+        e = obj_array_vectorize(
+                lambda el: discr.opposite_face_connection()(el),
+                i)
+
+    from grudge.symbolic.primitives import TracePair
+    return TracePair("int_faces", i, e)
+
 
 # vim: foldmethod=marker
diff --git a/grudge/execution.py b/grudge/execution.py
index a95b15e75e43917688f2c3766651e3018ff9d526..0104584ef28a8707f893b1c9d0ffd576c3453ade 100644
--- a/grudge/execution.py
+++ b/grudge/execution.py
@@ -22,12 +22,19 @@ OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
 THE SOFTWARE.
 """
 
-import six
+from typing import Optional, Union, Dict
+from numbers import Number
 import numpy as np
+
+from pytools import memoize_in
+from pytools.obj_array import make_obj_array
+
 import loopy as lp
 import pyopencl as cl
 import pyopencl.array  # noqa
-from pytools import memoize_in
+
+from meshmode.dof_array import DOFArray, thaw, flatten, unflatten
+from meshmode.array_context import ArrayContext, make_loopy_program
 
 import grudge.symbolic.mappers as mappers
 from grudge import sym
@@ -42,17 +49,20 @@ from loopy.version import LOOPY_USE_LANGUAGE_VERSION_2018_1  # noqa: F401
 MPI_TAG_SEND_TAGS = 1729
 
 
+ResultType = Union[DOFArray, Number]
+
+
 # {{{ exec mapper
 
 class ExecutionMapper(mappers.Evaluator,
         mappers.BoundOpMapperMixin,
         mappers.LocalOpReducerMixin):
-    def __init__(self, queue, context, bound_op):
+    def __init__(self, array_context, context, bound_op):
         super(ExecutionMapper, self).__init__(context)
         self.discrwb = bound_op.discrwb
         self.bound_op = bound_op
         self.function_registry = bound_op.function_registry
-        self.queue = queue
+        self.array_context = array_context
 
     # {{{ expression mappings
 
@@ -62,13 +72,14 @@ class ExecutionMapper(mappers.Evaluator,
 
         discr = self.discrwb.discr_from_dd(expr.dd)
 
-        result = discr.empty(self.queue, allocator=self.bound_op.allocator)
-        result.fill(1)
+        result = discr.empty(self.array_context)
+        for grp_ary in result:
+            grp_ary.fill(1)
         return result
 
     def map_node_coordinate_component(self, expr):
         discr = self.discrwb.discr_from_dd(expr.dd)
-        return discr.nodes()[expr.axis].with_queue(self.queue)
+        return thaw(self.array_context, discr.nodes()[expr.axis])
 
     def map_grudge_variable(self, expr):
         from numbers import Number
@@ -76,8 +87,9 @@ class ExecutionMapper(mappers.Evaluator,
         value = self.context[expr.name]
         if not expr.dd.is_scalar() and isinstance(value, Number):
             discr = self.discrwb.discr_from_dd(expr.dd)
-            ary = discr.empty(self.queue)
-            ary.fill(value)
+            ary = discr.empty(self.array_context)
+            for grp_ary in ary:
+                grp_ary.fill(value)
             value = ary
 
         return value
@@ -91,42 +103,41 @@ class ExecutionMapper(mappers.Evaluator,
             from numbers import Number
             if not dd.is_scalar() and isinstance(value, Number):
                 discr = self.discrwb.discr_from_dd(dd)
-                ary = discr.empty(self.queue)
-                ary.fill(value)
+                ary = discr.empty(self.array_context)
+                for grp_ary in ary:
+                    grp_ary.fill(value)
                 value = ary
         return value
 
     def map_call(self, expr):
         args = [self.rec(p) for p in expr.parameters]
-        return self.function_registry[expr.function.name](self.queue, *args)
+        return self.function_registry[expr.function.name](self.array_context, *args)
 
     # }}}
 
     # {{{ elementwise reductions
 
     def _map_elementwise_reduction(self, op_name, field_expr, dd):
-        @memoize_in(self, "elementwise_%s_knl" % op_name)
-        def knl():
-            knl = lp.make_kernel(
-                "{[el, idof, jdof]: 0<=el<nelements and 0<=idof, jdof<ndofs}",
+        @memoize_in(self, "elementwise_%s_prg" % op_name)
+        def prg():
+            return make_loopy_program(
+                "{[iel, idof, jdof]: 0<=iel<nelements and 0<=idof, jdof<ndofs}",
                 """
-                result[el, idof] = %s(jdof, operand[el, jdof])
+                result[iel, idof] = %s(jdof, operand[iel, jdof])
                 """ % op_name,
-                default_offset=lp.auto,
-                name="elementwise_%s_knl" % op_name)
-
-            return lp.tag_inames(knl, "el:g.0,idof:l.0")
+                name="grudge_elementwise_%s" % op_name)
 
         field = self.rec(field_expr)
         discr = self.discrwb.discr_from_dd(dd)
-        assert field.shape == (discr.nnodes,)
+        assert field.shape == (len(discr.groups),)
 
-        result = discr.empty(queue=self.queue, dtype=field.dtype,
-                allocator=self.bound_op.allocator)
+        result = discr.empty(self.array_context, dtype=field.entry_dtype)
         for grp in discr.groups:
-            knl()(self.queue,
-                    operand=grp.view(field),
-                    result=grp.view(result))
+            assert field[grp.index].shape == (grp.nelements, grp.nunit_dofs)
+            self.array_context.call_loopy(
+                    prg(),
+                    operand=field[grp.index],
+                    result=result[grp.index])
 
         return result
 
@@ -145,57 +156,101 @@ class ExecutionMapper(mappers.Evaluator,
 
     def map_nodal_sum(self, op, field_expr):
         # FIXME: Could allow array scalars
-        return cl.array.sum(self.rec(field_expr)).get()[()]
+        # FIXME: Fix CL-specific-ness
+        return sum([
+                cl.array.sum(grp_ary).get()[()]
+                for grp_ary in self.rec(field_expr)
+                ])
 
     def map_nodal_max(self, op, field_expr):
         # FIXME: Could allow array scalars
-        return cl.array.max(self.rec(field_expr)).get()[()]
+        # FIXME: Fix CL-specific-ness
+        return np.max([
+            cl.array.max(grp_ary).get()[()]
+            for grp_ary in self.rec(field_expr)])
 
     def map_nodal_min(self, op, field_expr):
         # FIXME: Could allow array scalars
-        return cl.array.min(self.rec(field_expr)).get()[()]
+        # FIXME: Fix CL-specific-ness
+        return np.min([
+            cl.array.min(grp_ary).get()[()]
+            for grp_ary in self.rec(field_expr)])
 
     # }}}
 
     def map_if(self, expr):
         bool_crit = self.rec(expr.condition)
 
+        if isinstance(bool_crit,  DOFArray):
+            # continues below
+            pass
+        elif isinstance(bool_crit,  np.number):
+            if bool_crit:
+                return self.rec(expr.then)
+            else:
+                return self.rec(expr.else_)
+        else:
+            raise TypeError(
+                "Expected criterion to be of type np.number or DOFArray")
+
+        assert isinstance(bool_crit,  DOFArray)
+        ngroups = len(bool_crit)
+
+        from pymbolic import var
+        iel = var("iel")
+        idof = var("idof")
+        subscript = (iel, idof)
+
         then = self.rec(expr.then)
         else_ = self.rec(expr.else_)
 
         import pymbolic.primitives as p
         var = p.Variable
 
-        i = var("i")
-        if isinstance(then,  pyopencl.array.Array):
-            sym_then = var("a")[i]
-        elif isinstance(then,  np.number):
+        if isinstance(then,  DOFArray):
+            sym_then = var("a")[subscript]
+
+            def get_then(igrp):
+                return then[igrp]
+        elif isinstance(then, np.number):
             sym_then = var("a")
+
+            def get_then(igrp):
+                return then
         else:
             raise TypeError(
-                "Expected parameter to be of type np.number or pyopencl.array.Array")
+                "Expected 'then' to be of type np.number or DOFArray")
 
-        if isinstance(else_,  pyopencl.array.Array):
-            sym_else = var("b")[i]
+        if isinstance(else_,  DOFArray):
+            sym_else = var("b")[subscript]
+
+            def get_else(igrp):
+                return else_[igrp]
         elif isinstance(else_,  np.number):
             sym_else = var("b")
+
+            def get_else(igrp):
+                return else_
         else:
             raise TypeError(
-                "Expected parameter to be of type np.number or pyopencl.array.Array")
+                "Expected 'else' to be of type np.number or DOFArray")
 
         @memoize_in(self.bound_op, "map_if_knl")
-        def knl():
-            knl = lp.make_kernel(
-                "{[i]: 0<=i<n}",
+        def knl(sym_then, sym_else):
+            return make_loopy_program(
+                "{[iel, idof]: 0<=iel<nelements and 0<=idof<nunit_dofs}",
                 [
-                    lp.Assignment(var("out")[i],
-                        p.If(var("crit")[i], sym_then, sym_else))
+                    lp.Assignment(var("out")[iel, idof],
+                        p.If(var("crit")[iel, idof], sym_then, sym_else))
                 ])
-            return lp.split_iname(knl, "i", 128, outer_tag="g.0", inner_tag="l.0")
-
-        evt, (out,) = knl()(self.queue, crit=bool_crit, a=then, b=else_)
 
-        return out
+        return DOFArray.from_list(self.array_context, [
+            self.array_context.call_loopy(
+                knl(sym_then, sym_else),
+                crit=bool_crit[igrp],
+                a=get_then(igrp),
+                b=get_else(igrp))
+            for igrp in range(ngroups)])
 
     # {{{ elementwise linear operators
 
@@ -211,25 +266,22 @@ class ExecutionMapper(mappers.Evaluator,
             return 0
 
         @memoize_in(self.bound_op, "elwise_linear_knl")
-        def knl():
-            knl = lp.make_kernel(
-                """{[k,i,j]:
-                    0<=k<nelements and
-                    0<=i<ndiscr_nodes_out and
+        def prg():
+            result = make_loopy_program(
+                """{[iel, idof, j]:
+                    0<=iel<nelements and
+                    0<=idof<ndiscr_nodes_out and
                     0<=j<ndiscr_nodes_in}""",
-                "result[k,i] = sum(j, mat[i, j] * vec[k, j])",
-                default_offset=lp.auto, name="diff")
+                "result[iel, idof] = sum(j, mat[idof, j] * vec[iel, j])",
+                name="diff")
 
-            knl = lp.split_iname(knl, "i", 16, inner_tag="l.0")
-            knl = lp.tag_array_axes(knl, "mat", "stride:auto,stride:auto")
-            return lp.tag_inames(knl, dict(k="g.0"))
+            result = lp.tag_array_axes(result, "mat", "stride:auto,stride:auto")
+            return result
 
         in_discr = self.discrwb.discr_from_dd(op.dd_in)
         out_discr = self.discrwb.discr_from_dd(op.dd_out)
 
-        result = out_discr.empty(
-                queue=self.queue,
-                dtype=field.dtype, allocator=self.bound_op.allocator)
+        result = out_discr.empty(self.array_context, dtype=field.entry_dtype)
 
         for in_grp, out_grp in zip(in_discr.groups, out_discr.groups):
 
@@ -237,32 +289,34 @@ class ExecutionMapper(mappers.Evaluator,
             try:
                 matrix = self.bound_op.operator_data_cache[cache_key]
             except KeyError:
-                matrix = (
-                    cl.array.to_device(
-                        self.queue,
-                        np.asarray(op.matrix(out_grp, in_grp), dtype=field.dtype))
-                    .with_queue(None))
+                matrix = self.array_context.freeze(
+                        self.array_context.from_numpy(
+                            np.asarray(
+                                op.matrix(out_grp, in_grp),
+                                dtype=field.entry_dtype)))
 
                 self.bound_op.operator_data_cache[cache_key] = matrix
 
-            knl()(self.queue, mat=matrix, result=out_grp.view(result),
-                    vec=in_grp.view(field))
+            self.array_context.call_loopy(
+                    prg(),
+                    mat=matrix,
+                    result=result[out_grp.index],
+                    vec=field[in_grp.index])
 
         return result
 
     def map_interpolation(self, op, field_expr):
         conn = self.discrwb.connection_from_dds(op.dd_in, op.dd_out)
-        return conn(self.queue, self.rec(field_expr)).with_queue(self.queue)
+        return conn(self.rec(field_expr))
 
     def map_opposite_partition_face_swap(self, op, field_expr):
         assert op.dd_in == op.dd_out
         bdry_conn = self.discrwb.get_distributed_boundary_swap_connection(op.dd_in)
         remote_bdry_vec = self.rec(field_expr)  # swapped by RankDataSwapAssign
-        return bdry_conn(self.queue, remote_bdry_vec).with_queue(self.queue)
+        return bdry_conn(remote_bdry_vec)
 
     def map_opposite_interior_face_swap(self, op, field_expr):
-        return self.discrwb.opposite_face_connection()(
-                self.queue, self.rec(field_expr)).with_queue(self.queue)
+        return self.discrwb.opposite_face_connection()(self.rec(field_expr))
 
     # }}}
 
@@ -276,26 +330,23 @@ class ExecutionMapper(mappers.Evaluator,
             return 0
 
         @memoize_in(self.bound_op, "face_mass_knl")
-        def knl():
-            knl = lp.make_kernel(
-                """{[k,i,f,j]:
-                    0<=k<nelements and
+        def prg():
+            return make_loopy_program(
+                """{[iel,idof,f,j]:
+                    0<=iel<nelements and
                     0<=f<nfaces and
-                    0<=i<nvol_nodes and
+                    0<=idof<nvol_nodes and
                     0<=j<nface_nodes}""",
-                "result[k,i] = sum(f, sum(j, mat[i, f, j] * vec[f, k, j]))",
-                default_offset=lp.auto, name="face_mass")
-
-            knl = lp.split_iname(knl, "i", 16, inner_tag="l.0")
-            return lp.tag_inames(knl, dict(k="g.0"))
+                """
+                result[iel,idof] = sum(f, sum(j, mat[idof, f, j] * vec[f, iel, j]))
+                """,
+                name="face_mass")
 
         all_faces_conn = self.discrwb.connection_from_dds("vol", op.dd_in)
         all_faces_discr = all_faces_conn.to_discr
         vol_discr = all_faces_conn.from_discr
 
-        result = vol_discr.empty(
-                queue=self.queue,
-                dtype=field.dtype, allocator=self.bound_op.allocator)
+        result = vol_discr.empty(self.array_context, dtype=field.entry_dtype)
 
         assert len(all_faces_discr.groups) == len(vol_discr.groups)
 
@@ -307,16 +358,18 @@ class ExecutionMapper(mappers.Evaluator,
             try:
                 matrix = self.bound_op.operator_data_cache[cache_key]
             except KeyError:
-                matrix = op.matrix(afgrp, volgrp, field.dtype)
-                matrix = (
-                        cl.array.to_device(self.queue, matrix)
-                        .with_queue(None))
+                matrix = op.matrix(afgrp, volgrp, field.entry_dtype)
+                matrix = self.array_context.freeze(
+                        self.array_context.from_numpy(matrix))
 
                 self.bound_op.operator_data_cache[cache_key] = matrix
 
-            input_view = afgrp.view(field).reshape(
-                    nfaces, volgrp.nelements, afgrp.nunit_nodes)
-            knl()(self.queue, mat=matrix, result=volgrp.view(result),
+            input_view = field[afgrp.index].reshape(
+                    nfaces, volgrp.nelements, afgrp.nunit_dofs)
+            self.array_context.call_loopy(
+                    prg(),
+                    mat=matrix,
+                    result=result[volgrp.index],
                     vec=input_view)
 
         return result
@@ -332,15 +385,16 @@ class ExecutionMapper(mappers.Evaluator,
                 sym.DD_VOLUME,
                 sym.DOFDesc(expr.dd.domain_tag))
 
-        field = face_discr.empty(self.queue,
-                dtype=self.discrwb.real_dtype,
-                allocator=self.bound_op.allocator)
-        field.fill(1)
+        field = face_discr.empty(self.array_context, dtype=self.discrwb.real_dtype)
+        for grp_ary in field:
+            grp_ary.fill(1)
 
-        for grp in all_faces_conn.groups:
+        for igrp, grp in enumerate(all_faces_conn.groups):
             for batch in grp.batches:
-                i = batch.to_element_indices.with_queue(self.queue)
-                field[i] = (2.0 * (batch.to_element_face % 2) - 1.0) * field[i]
+                i = self.array_context.thaw(batch.to_element_indices)
+                grp_field = field[igrp].reshape(-1)
+                grp_field[i] = \
+                        (2.0 * (batch.to_element_face % 2) - 1.0) * grp_field[i]
 
         return field
 
@@ -349,7 +403,7 @@ class ExecutionMapper(mappers.Evaluator,
     # {{{ instruction execution functions
 
     def map_insn_rank_data_swap(self, insn, profile_data=None):
-        local_data = self.rec(insn.field).get(self.queue)
+        local_data = self.array_context.to_numpy(flatten(self.rec(insn.field)))
         comm = self.discrwb.mpi_communicator
 
         # print("Sending data to rank %d with tag %d"
@@ -360,31 +414,60 @@ class ExecutionMapper(mappers.Evaluator,
         recv_req = comm.Irecv(remote_data_host, insn.i_remote_rank, insn.recv_tag)
 
         return [], [
-                MPIRecvFuture(recv_req, insn.name, remote_data_host, self.queue),
+                MPIRecvFuture(
+                    array_context=self.array_context,
+                    bdry_discr=self.discrwb.discr_from_dd(insn.dd_out),
+                    recv_req=recv_req,
+                    insn_name=insn.name,
+                    remote_data_host=remote_data_host),
                 MPISendFuture(send_req)]
 
     def map_insn_loopy_kernel(self, insn, profile_data=None):
-        kwargs = {}
         kdescr = insn.kernel_descriptor
-        for name, expr in six.iteritems(kdescr.input_mappings):
-            kwargs[name] = self.rec(expr)
-
         discr = self.discrwb.discr_from_dd(kdescr.governing_dd)
+
+        dof_array_kwargs = {}
+        other_kwargs = {}
+
+        for name, expr in kdescr.input_mappings.items():
+            v = self.rec(expr)
+            if isinstance(v, DOFArray):
+                dof_array_kwargs[name] = v
+            else:
+                other_kwargs[name] = v
+
         for name in kdescr.scalar_args():
-            v = kwargs[name]
+            v = other_kwargs[name]
             if isinstance(v, (int, float)):
-                kwargs[name] = discr.real_dtype.type(v)
+                other_kwargs[name] = discr.real_dtype.type(v)
             elif isinstance(v, complex):
-                kwargs[name] = discr.complex_dtype.type(v)
+                other_kwargs[name] = discr.complex_dtype.type(v)
             elif isinstance(v, np.number):
                 pass
             else:
                 raise ValueError("unrecognized scalar type for variable '%s': %s"
                         % (name, type(v)))
 
-        kwargs["grdg_n"] = discr.nnodes
-        evt, result_dict = kdescr.loopy_kernel(self.queue, **kwargs)
-        return list(result_dict.items()), []
+        result = {}
+        for grp in discr.groups:
+            kwargs = other_kwargs.copy()
+            kwargs["nelements"] = grp.nelements
+            kwargs["nunit_dofs"] = grp.nunit_dofs
+
+            for name, ary in dof_array_kwargs.items():
+                kwargs[name] = ary[grp.index]
+
+            knl_result = self.array_context.call_loopy(
+                    kdescr.loopy_kernel, **kwargs)
+
+            for name, val in knl_result.items():
+                result.setdefault(name, []).append(val)
+
+        result = {
+                name: DOFArray.from_list(self.array_context, val)
+                for name, val in result.items()}
+
+        return list(result.items()), []
 
     def map_insn_assign(self, insn, profile_data=None):
         return [(name, self.rec(expr))
@@ -413,35 +496,35 @@ class ExecutionMapper(mappers.Evaluator,
 
         assert repr_op.dd_in.domain_tag == repr_op.dd_out.domain_tag
 
-        @memoize_in(self.discrwb, "reference_derivative_knl")
-        def knl():
-            knl = lp.make_kernel(
-                """{[imatrix,k,i,j]:
+        @memoize_in(self.discrwb, "reference_derivative_prg")
+        def prg(nmatrices):
+            result = make_loopy_program(
+                """{[imatrix, iel, idof, j]:
                     0<=imatrix<nmatrices and
-                    0<=k<nelements and
-                    0<=i<nunit_nodes_out and
+                    0<=iel<nelements and
+                    0<=idof<nunit_nodes_out and
                     0<=j<nunit_nodes_in}""",
                 """
-                result[imatrix, k, i] = sum(
-                        j, diff_mat[imatrix, i, j] * vec[k, j])
+                result[imatrix, iel, idof] = sum(
+                        j, diff_mat[imatrix, idof, j] * vec[iel, j])
                 """,
-                default_offset=lp.auto, name="diff")
+                name="diff")
 
-            knl = lp.split_iname(knl, "i", 16, inner_tag="l.0")
-            return lp.tag_inames(knl, dict(k="g.0"))
+            result = lp.fix_parameters(result, nmatrices=nmatrices)
+            result = lp.tag_inames(result, "imatrix: unr")
+            result = lp.tag_array_axes(result, "result", "sep,c,c")
+            return result
 
         noperators = len(insn.operators)
 
         in_discr = self.discrwb.discr_from_dd(repr_op.dd_in)
         out_discr = self.discrwb.discr_from_dd(repr_op.dd_out)
 
-        result = out_discr.empty(
-                queue=self.queue,
-                dtype=field.dtype, extra_dims=(noperators,),
-                allocator=self.bound_op.allocator)
+        result = make_obj_array([
+            out_discr.empty(self.array_context, dtype=field.entry_dtype)
+            for idim in range(noperators)])
 
         for in_grp, out_grp in zip(in_discr.groups, out_discr.groups):
-
             if in_grp.nelements == 0:
                 continue
 
@@ -449,13 +532,18 @@ class ExecutionMapper(mappers.Evaluator,
 
             # FIXME: Should transfer matrices to device and cache them
             matrices_ary = np.empty((
-                noperators, out_grp.nunit_nodes, in_grp.nunit_nodes))
+                noperators, out_grp.nunit_dofs, in_grp.nunit_dofs))
             for i, op in enumerate(insn.operators):
                 matrices_ary[i] = matrices[op.rst_axis]
+            matrices_ary_dev = self.array_context.from_numpy(matrices_ary)
 
-            knl()(self.queue,
-                    diff_mat=matrices_ary,
-                    result=out_grp.view(result), vec=in_grp.view(field))
+            self.array_context.call_loopy(
+                    prg(noperators),
+                    diff_mat=matrices_ary_dev,
+                    result=make_obj_array([
+                        result[iop][out_grp.index]
+                        for iop in range(noperators)
+                        ]), vec=field[in_grp.index])
 
         return [(name, result[i]) for i, name in enumerate(insn.names)], []
 
@@ -467,18 +555,22 @@ class ExecutionMapper(mappers.Evaluator,
 # {{{ futures
 
 class MPIRecvFuture(object):
-    def __init__(self, recv_req, insn_name, remote_data_host, queue):
+    def __init__(self, array_context, bdry_discr, recv_req, insn_name,
+            remote_data_host):
+        self.array_context = array_context
+        self.bdry_discr = bdry_discr
         self.receive_request = recv_req
         self.insn_name = insn_name
         self.remote_data_host = remote_data_host
-        self.queue = queue
 
     def is_ready(self):
         return self.receive_request.Test()
 
     def __call__(self):
         self.receive_request.Wait()
-        remote_data = cl.array.to_device(self.queue, self.remote_data_host)
+        actx = self.array_context
+        remote_data = unflatten(self.array_context, self.bdry_discr,
+                actx.from_numpy(self.remote_data_host))
         return [(self.insn_name, remote_data)], []
 
 
@@ -499,16 +591,14 @@ class MPISendFuture(object):
 # {{{ bound operator
 
 class BoundOperator(object):
-
     def __init__(self, discrwb, discr_code, eval_code, debug_flags,
-            function_registry, exec_mapper_factory, allocator=None):
+            function_registry, exec_mapper_factory):
         self.discrwb = discrwb
         self.discr_code = discr_code
         self.eval_code = eval_code
         self.operator_data_cache = {}
         self.debug_flags = debug_flags
         self.function_registry = function_registry
-        self.allocator = allocator
         self.exec_mapper_factory = exec_mapper_factory
 
     def __str__(self):
@@ -523,16 +613,50 @@ class BoundOperator(object):
                 + sep
                 + str(self.eval_code))
 
-    def __call__(self, queue, profile_data=None, log_quantities=None, **context):
-        import pyopencl.array as cl_array
-
-        def replace_queue(a):
-            if isinstance(a, cl_array.Array):
-                return a.with_queue(queue)
+    def __call__(self, array_context: Optional[ArrayContext] = None,
+            *, profile_data=None, log_quantities=None, **context):
+        """
+        :arg array_context: only needs to be supplied if no instances of
+            :class:`~meshmode.dof_array.DOFArray` with a
+            :class:`~meshmode.array_context.ArrayContext`
+            are supplied as part of *context*.
+        """
+
+        # {{{ figure array context
+
+        array_contexts = []
+        if array_context is not None:
+            if not isinstance(array_context, ArrayContext):
+                raise TypeError(
+                        "first positional argument (if supplied) must be "
+                        "an ArrayContext")
+
+            array_contexts.append(array_context)
+        del array_context
+
+        def look_for_array_contexts(ary):
+            if isinstance(ary, DOFArray):
+                if ary.array_context is not None:
+                    array_contexts.append(ary.array_context)
+            elif isinstance(ary, np.ndarray) and ary.dtype.char == "O":
+                for idx in np.ndindex(ary.shape):
+                    look_for_array_contexts(ary[idx])
             else:
-                return a
+                pass
+
+        for key, val in context.items():
+            look_for_array_contexts(val)
 
-        from pytools.obj_array import with_object_array_or_scalar
+        if array_contexts:
+            from pytools import is_single_valued
+            if not is_single_valued(array_contexts):
+                raise ValueError("arguments do not agree on an array context")
+
+            array_context = array_contexts[0]
+        else:
+            raise ValueError("no array context given or available from arguments")
+
+        # }}}
 
         # {{{ discrwb-scope evaluation
 
@@ -541,18 +665,15 @@ class BoundOperator(object):
                     self.discrwb._discr_scoped_subexpr_name_to_value)
                 for result_var in self.discr_code.result):
             # need to do discrwb-scope evaluation
-            discrwb_eval_context = {}
+            discrwb_eval_context: Dict[str, ResultType] = {}
             self.discr_code.execute(
-                    self.exec_mapper_factory(queue, discrwb_eval_context, self))
+                    self.exec_mapper_factory(
+                        array_context, discrwb_eval_context, self))
 
         # }}}
 
-        new_context = {}
-        for name, var in six.iteritems(context):
-            new_context[name] = with_object_array_or_scalar(replace_queue, var)
-
         return self.eval_code.execute(
-                self.exec_mapper_factory(queue, new_context, self),
+                self.exec_mapper_factory(array_context, context, self),
                 profile_data=profile_data,
                 log_quantities=log_quantities)
 
@@ -657,7 +778,7 @@ def process_sym_operator(discrwb, sym_operator, post_bind_mapper=None,
 def bind(discr, sym_operator, post_bind_mapper=lambda x: x,
         function_registry=base_function_registry,
         exec_mapper_factory=ExecutionMapper,
-        debug_flags=frozenset(), allocator=None):
+        debug_flags=frozenset()):
     # from grudge.symbolic.mappers import QuadratureUpsamplerRemover
     # sym_operator = QuadratureUpsamplerRemover(self.quad_min_degrees)(
     #         sym_operator)
@@ -685,8 +806,7 @@ def bind(discr, sym_operator, post_bind_mapper=lambda x: x,
     bound_op = BoundOperator(discr, discr_code, eval_code,
             function_registry=function_registry,
             exec_mapper_factory=exec_mapper_factory,
-            debug_flags=debug_flags,
-            allocator=allocator)
+            debug_flags=debug_flags)
 
     if "dump_op_code" in debug_flags:
         from pytools.debug import open_unique_debug_file
diff --git a/grudge/function_registry.py b/grudge/function_registry.py
index 4d521560d8e6c83a8e9fa238b3c46955dc62e7a7..14970c8a6aed0cd2ea63f71ba4337accb33b81cd 100644
--- a/grudge/function_registry.py
+++ b/grudge/function_registry.py
@@ -26,11 +26,10 @@ THE SOFTWARE.
 """
 
 
-import loopy as lp
 import numpy as np
 
 from loopy.version import LOOPY_USE_LANGUAGE_VERSION_2018_1  # noqa
-from pytools import RecordWithoutPickling, memoize_in
+from pytools import RecordWithoutPickling
 
 
 # {{{ helpers
@@ -82,76 +81,29 @@ class Function(RecordWithoutPickling):
         raise NotImplementedError
 
 
-class CElementwiseUnaryFunction(Function):
-
+class CElementwiseFunction(Function):
     supports_codegen = True
 
+    def __init__(self, identifier, nargs):
+        super().__init__(identifier=identifier, nargs=nargs)
+
     def get_result_dofdesc(self, arg_dds):
-        assert len(arg_dds) == 1
+        assert len(arg_dds) == self.nargs
         return arg_dds[0]
 
-    def __call__(self, queue, arg):
+    def __call__(self, array_context, *args):
         func_name = self.identifier
-        if should_use_numpy(arg):
+        from pytools import single_valued
+        if single_valued(should_use_numpy(arg) for arg in args):
             func = getattr(np, func_name)
-            return func(arg)
+            return func(*args)
 
-        cached_name = "map_call_knl_"
-
-        from pymbolic.primitives import Variable
-        i = Variable("i")
-
-        if self.identifier == "fabs":  # FIXME
+        if func_name == "fabs":  # FIXME
             # Loopy has a type-adaptive "abs", but no "fabs".
             func_name = "abs"
 
-        cached_name += func_name
-
-        @memoize_in(self, cached_name)
-        def knl():
-            knl = lp.make_kernel(
-                "{[i]: 0<=i<n}",
-                [
-                    lp.Assignment(Variable("out")[i],
-                        Variable(func_name)(Variable("a")[i]))
-                ], default_offset=lp.auto)
-            return lp.split_iname(knl, "i", 128, outer_tag="g.0", inner_tag="l.0")
-
-        evt, (out,) = knl()(queue, a=arg)
-        return out
-
-
-class CElementwiseBinaryFunction(Function):
-    supports_codegen = True
-
-    def get_result_dofdesc(self, arg_dds):
-        assert len(arg_dds) == 2
-        from pytools import single_valued
-        return single_valued(arg_dds)
-
-    def __call__(self, queue, arg0, arg1):
-        func_name = self.identifier
-        if should_use_numpy(arg0) and should_use_numpy(arg1):
-            func = getattr(np, cl_to_numpy_function_name(func_name))
-            return func(arg0, arg1)
-
-        from pymbolic.primitives import Variable
-
-        @memoize_in(self, "map_call_knl_%s" % func_name)
-        def knl():
-            i = Variable("i")
-            knl = lp.make_kernel(
-                "{[i]: 0<=i<n}",
-                [
-                    lp.Assignment(
-                        Variable("out")[i],
-                        Variable(func_name)(Variable("a")[i], Variable("b")[i]))
-                ], default_offset=lp.auto)
-
-            return lp.split_iname(knl, "i", 128, outer_tag="g.0", inner_tag="l.0")
-
-        evt, (out,) = knl()(queue, a=arg0, b=arg1)
-        return out
+        sfunc = array_context.special_func(func_name)
+        return sfunc(*args)
 
 
 class CBesselFunction(Function):
@@ -173,8 +125,8 @@ class FixedDOFDescExternalFunction(Function):
                 implementation=implementation,
                 dd=dd)
 
-    def __call__(self, queue, *args, **kwargs):
-        return self.implementation(queue, *args, **kwargs)
+    def __call__(self, array_context, *args, **kwargs):
+        return self.implementation(array_context, *args, **kwargs)
 
     def get_result_dofdesc(self, arg_dds):
         return self.dd
@@ -220,12 +172,12 @@ class FunctionRegistry(RecordWithoutPickling):
 def _make_bfr():
     bfr = FunctionRegistry()
 
-    bfr = bfr.register(CElementwiseUnaryFunction("sqrt"))
-    bfr = bfr.register(CElementwiseUnaryFunction("exp"))
-    bfr = bfr.register(CElementwiseUnaryFunction("fabs"))
-    bfr = bfr.register(CElementwiseUnaryFunction("sin"))
-    bfr = bfr.register(CElementwiseUnaryFunction("cos"))
-    bfr = bfr.register(CElementwiseBinaryFunction("atan2"))
+    bfr = bfr.register(CElementwiseFunction("sqrt", 1))
+    bfr = bfr.register(CElementwiseFunction("exp", 1))
+    bfr = bfr.register(CElementwiseFunction("fabs", 1))
+    bfr = bfr.register(CElementwiseFunction("sin", 1))
+    bfr = bfr.register(CElementwiseFunction("cos", 1))
+    bfr = bfr.register(CElementwiseFunction("atan2", 2))
     bfr = bfr.register(CBesselFunction("bessel_j"))
     bfr = bfr.register(CBesselFunction("bessel_y"))
 
diff --git a/grudge/shortcuts.py b/grudge/shortcuts.py
index 475f72d24be204627b114813493ce623c431f9c4..94b8fb2ac6674667f8753aed501fbcc1806a1ffc 100644
--- a/grudge/shortcuts.py
+++ b/grudge/shortcuts.py
@@ -25,9 +25,6 @@ THE SOFTWARE.
 """
 
 
-import pyopencl as cl
-
-
 def set_up_rk4(field_var_name, dt, fields, rhs, t_start=0):
     from leap.rk import LSRK4MethodBuilder
     from dagrt.codegen import PythonCodeGenerator
@@ -46,14 +43,14 @@ def set_up_rk4(field_var_name, dt, fields, rhs, t_start=0):
 
 def make_visualizer(discrwb, vis_order):
     from meshmode.discretization.visualization import make_visualizer
-    with cl.CommandQueue(discrwb.cl_context) as queue:
-        return make_visualizer(queue, discrwb.discr_from_dd("vol"), vis_order)
+    return make_visualizer(
+            discrwb._setup_actx,
+            discrwb.discr_from_dd("vol"), vis_order)
 
 
 def make_boundary_visualizer(discrwb, vis_order):
     from meshmode.discretization.visualization import make_visualizer
     from grudge import sym
-    with cl.CommandQueue(discrwb.cl_context) as queue:
-        return make_visualizer(
-                queue, discrwb.discr_from_dd(sym.BTAG_ALL),
-                vis_order)
+    return make_visualizer(
+            discrwb._setup_actx, discrwb.discr_from_dd(sym.BTAG_ALL),
+            vis_order)
diff --git a/grudge/symbolic/compiler.py b/grudge/symbolic/compiler.py
index 51eb5202acaa1d1ffa6c4cbc8b293eef4a25e54e..cedd152e653418e6776ab36421050bf880368fa4 100644
--- a/grudge/symbolic/compiler.py
+++ b/grudge/symbolic/compiler.py
@@ -87,7 +87,7 @@ class LoopyKernelDescriptor(object):
         import loopy as lp
         return [arg.name for arg in self.loopy_kernel.args
                 if isinstance(arg, lp.ValueArg)
-                and arg.name != "grdg_n"]
+                and arg.name not in ["nelements", "nunit_dofs"]]
 
 
 class LoopyKernelInstruction(Instruction):
@@ -845,13 +845,11 @@ def is_function_loopyable(function, function_registry):
 
 
 class ToLoopyExpressionMapper(mappers.IdentityMapper):
-    def __init__(self, dd_inference_mapper, temp_names, iname):
+    def __init__(self, dd_inference_mapper, temp_names, subscript):
         self.dd_inference_mapper = dd_inference_mapper
         self.function_registry = dd_inference_mapper.function_registry
         self.temp_names = temp_names
-        self.iname = iname
-        from pymbolic import var
-        self.iname_expr = var(iname)
+        self.subscript = subscript
 
         self.expr_to_name = {}
         self.used_names = set()
@@ -887,7 +885,7 @@ class ToLoopyExpressionMapper(mappers.IdentityMapper):
             return var(name)
         else:
             self.non_scalar_vars.append(name)
-            return var(name)[self.iname_expr]
+            return var(name)[self.subscript]
 
     def map_variable(self, expr):
         return self.map_variable_ref_expr(expr, expr.name)
@@ -1004,16 +1002,19 @@ class ToLoopyInstructionMapper(object):
                         insn.exprs[0], self.function_registry))):
             return insn
 
-        iname = "grdg_i"
-        size_name = "grdg_n"
+        # FIXME: These names and the size names could clash with user-given names.
+        # Need better metadata tracking in loopy.
+        iel = "iel"
+        idof = "idof"
 
         temp_names = [
                 name
                 for name, dnr in zip(insn.names, insn.do_not_return)
                 if dnr]
 
+        from pymbolic import var
         expr_mapper = ToLoopyExpressionMapper(
-                self.dd_inference_mapper, temp_names, iname)
+                self.dd_inference_mapper, temp_names, (var(iel), var(idof)))
         insns = []
 
         import loopy as lp
@@ -1033,17 +1034,21 @@ class ToLoopyInstructionMapper(object):
             return insn
 
         knl = lp.make_kernel(
-                "{[%s]: 0 <= %s < %s}" % (iname, iname, size_name),
+                "{[%(iel)s, %(idof)s]: "
+                "0 <= %(iel)s < nelements and 0 <= %(idof)s < nunit_dofs}"
+                % {"iel": iel, "idof": idof},
                 insns,
-                default_offset=lp.auto,
 
                 name="grudge_assign_%d" % self.insn_count,
+
                 # Single-insn kernels may have their no_sync_with resolve to an
                 # empty set, that's OK.
-                options=lp.Options(check_dep_resolution=False))
-
-        knl = lp.set_options(knl, return_dict=True)
-        knl = lp.split_iname(knl, iname, 128, outer_tag="g.0", inner_tag="l.0")
+                options=lp.Options(
+                    check_dep_resolution=False,
+                    return_dict=True,
+                    no_numpy=True,
+                    )
+                )
 
         self.insn_count += 1