diff --git a/examples/advection/gas_dynamics_initials.py b/examples/advection/gas_dynamics_initials.py
new file mode 100644
index 0000000000000000000000000000000000000000..6f4e4f4164c4484cd0b1ded50fb798a750cda72e
--- /dev/null
+++ b/examples/advection/gas_dynamics_initials.py
@@ -0,0 +1,172 @@
+# Copyright (C) 2008 Andreas Kloeckner
+#
+# This program is free software: you can redistribute it and/or modify
+# it under the terms of the GNU General Public License as published by
+# the Free Software Foundation, either version 3 of the License, or
+# (at your option) any later version.
+#
+# This program is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+# GNU General Public License for more details.
+#
+# You should have received a copy of the GNU General Public License
+# along with this program.  If not, see <http://www.gnu.org/licenses/>.
+
+
+
+
+from __future__ import division
+from __future__ import absolute_import
+import numpy
+import numpy.linalg as la
+from six.moves import range
+from grudge import sym
+
+
+
+
+class UniformMachFlow:
+    def __init__(self, mach=0.1, p=1, rho=1, reynolds=100,
+            gamma=1.4, prandtl=0.72, char_length=1, spec_gas_const=287.1,
+            angle_of_attack=None, direction=None, gaussian_pulse_at=None,
+            pulse_magnitude=0.1):
+        """
+        :param direction: is a vector indicating the direction of the
+          flow. Only one of angle_of_attack and direction may be
+          specified. Only the direction, not the magnitude, of
+          direction is taken into account.
+
+        :param angle_of_attack: if not None, specifies the angle of
+          the flow along the Y axis, where the flow is
+          directed along the X axis.
+        """
+        if angle_of_attack is not None and direction is not None:
+            raise ValueError("Only one of angle_of_attack and "
+                    "direction may be specified.")
+
+        if angle_of_attack is None and direction is None:
+            angle_of_attack = 0
+
+        if direction is not None:
+            self.direction = direction/la.norm(direction)
+        else:
+            self.direction = None
+
+        self.mach = mach
+        self.p = p
+        self.rho = rho
+
+        self.gamma = gamma
+        self.prandtl = prandtl
+        self.reynolds = reynolds
+        self.length = char_length
+        self.spec_gas_const = spec_gas_const
+
+        self.angle_of_attack = angle_of_attack
+
+        self.gaussian_pulse_at = gaussian_pulse_at
+        self.pulse_magnitude = pulse_magnitude
+
+        self.c = (self.gamma * p / rho)**0.5
+        u = self.velocity = mach * self.c
+        self.e = p / (self.gamma - 1) + rho / 2 * u**2
+
+        if numpy.isinf(self.reynolds):
+            self.mu = 0
+        else:
+            self.mu = u * self.length * rho / self.reynolds
+
+    def direction_vector(self, dimensions):
+        # this must be done here because dimensions is not known above
+        if self.direction is None:
+            assert self.angle_of_attack is not None
+            direction = numpy.zeros(dimensions, dtype=numpy.float64)
+            direction[0] = numpy.cos(
+                    self.angle_of_attack / 180. * numpy.pi)
+            direction[1] = numpy.sin(
+                    self.angle_of_attack / 180. * numpy.pi)
+            return direction
+        else:
+            return self.direction
+
+    def __call__(self, t, x_vec):
+        ones = numpy.ones_like(x_vec[0])
+        rho_field = ones*self.rho
+
+        if self.gaussian_pulse_at is not None:
+            rel_to_pulse = [x_vec[i] - self.gaussian_pulse_at[i]
+                    for i in range(len(x_vec))]
+            rho_field +=  self.pulse_magnitude * self.rho * numpy.exp(
+                - sum(rtp_i**2 for rtp_i in rel_to_pulse)/2)
+
+        direction = self.direction_vector(x_vec.shape[0])
+
+        from grudge.tools import make_obj_array
+        u_field = make_obj_array([ones*self.velocity*dir_i
+            for dir_i in direction])
+
+        from grudge.tools import join_fields
+        return join_fields(rho_field, self.e*ones, self.rho*u_field)
+
+    def volume_interpolant(self, t, discr):
+        return discr.convert_volume(
+                        self(t, discr.nodes.T),
+                        kind=discr.compute_kind,
+                        dtype=discr.default_scalar_type)
+
+    def boundary_interpolant(self, t, discr, tag):
+        return discr.convert_boundary(
+                        self(t, discr.get_boundary(tag).nodes.T),
+                         tag=tag, kind=discr.compute_kind,
+                         dtype=discr.default_scalar_type)
+
+
+
+
+
+
+
+class Vortex:
+    def __init__(self):
+        self.beta = 5
+        self.gamma = 1.4
+        self.center = numpy.array([5, 0])
+        self.velocity = numpy.array([0.3, 0.3])
+        self.final_time = 0.5
+
+        self.mu = 0
+        self.prandtl = 0.72
+        self.spec_gas_const = 287.1
+
+    def __call__(self):
+        x_vec = sym.nodes(2)
+        t = sym.var("t", sym.DD_SCALAR)
+        vortex_loc = self.center + t*self.velocity
+
+        # coordinates relative to vortex center
+        x_rel = x_vec[0] - vortex_loc[0]
+        y_rel = x_vec[1] - vortex_loc[1]
+
+        # Y.C. Zhou, G.W. Wei / Journal of Computational Physics 189 (2003) 159
+        # also JSH/TW Nodal DG Methods, p. 209
+
+        from math import pi
+        r = sym.sqrt(x_rel**2+y_rel**2)
+        expterm = self.beta*sym.exp(1-r**2)
+        u = self.velocity[0] - expterm*y_rel/(2*pi)
+        v = self.velocity[1] + expterm*x_rel/(2*pi)
+        rho = (1-(self.gamma-1)/(16*self.gamma*pi**2)*expterm**2)**(1/(self.gamma-1))
+        p = rho**self.gamma
+
+        e = p/(self.gamma-1) + rho/2*(u**2+v**2)
+
+        from pytools.obj_array import join_fields
+        return join_fields(rho, e, rho*u, rho*v)
+
+    def volume_interpolant(self):
+        return self()
+
+
+
+
diff --git a/examples/advection/new.py b/examples/advection/new.py
new file mode 100644
index 0000000000000000000000000000000000000000..6a8566e064120009c3dff0028e7114e3b7219d5c
--- /dev/null
+++ b/examples/advection/new.py
@@ -0,0 +1,138 @@
+from __future__ import division, absolute_import
+
+__copyright__ = "Copyright (C) 2017 Bogdan Enache"
+
+__license__ = """
+Permission is hereby granted, free of charge, to any person obtaining a copy
+of this software and associated documentation files (the "Software"), to deal
+in the Software without restriction, including without limitation the rights
+to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+copies of the Software, and to permit persons to whom the Software is
+furnished to do so, subject to the following conditions:
+
+The above copyright notice and this permission notice shall be included in
+all copies or substantial portions of the Software.
+
+THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
+THE SOFTWARE.
+"""
+
+
+import numpy as np
+import pyopencl as cl  # noqa
+import pyopencl.array  # noqa
+import pyopencl.clmath  # noqa
+
+import pytest  # noqa
+
+from pyopencl.tools import (  # noqa
+                pytest_generate_tests_for_pyopencl as pytest_generate_tests)
+
+import logging
+logger = logging.getLogger(__name__)
+
+from grudge import sym, bind, DGDiscretizationWithBoundaries
+from pytools.obj_array import join_fields
+
+
+def main(write_output=True, order=4):
+    cl_ctx = cl.create_some_context()
+    queue = cl.CommandQueue(cl_ctx)
+
+    dim = 2
+
+    resolution = 15
+    from meshmode.mesh.generation import generate_regular_rect_mesh
+    mesh = generate_regular_rect_mesh(a=(0, -5), b=(10, 5),
+            n=(resolution, resolution), order=order)
+
+    h = 1/resolution
+
+
+    source_center = np.array([3, -2])
+    source_center2 = np.array([6, 2])
+    source_width = 0.5
+
+    sym_x = sym.nodes(2)
+    sym_source_center_dist = sym_x - source_center
+    sym_source_center_dist2 = sym_x - source_center2
+
+
+    def f_gaussian(x):
+        return join_fields(sym.exp(
+                    -np.dot(sym_source_center_dist, sym_source_center_dist)
+                    / source_width**2) ,sym.exp(
+                    -np.dot(sym_source_center_dist2, sym_source_center_dist2)
+                    / source_width**2) )
+
+    from gas_dynamics_initials import Vortex
+    flow = Vortex()
+    bund = sym.interp("vol", sym.BTAG_ALL)(flow.volume_interpolant()[0:2])
+
+    
+
+    discr = DGDiscretizationWithBoundaries(cl_ctx, mesh, order=order)
+
+    def sym_operator():
+        u = sym.make_sym_array("u", 2)
+        def mflux(pair):
+            normal = sym.normal(pair.dd, dim)
+            return pair.avg*np.dot(normal, [1,1])
+
+        def flux(pair):
+            return sym.interp(pair.dd, "all_faces")(
+                    mflux(pair))
+        stiff_t = sym.stiffness_t(dim)
+        return sym.InverseMassOperator()(
+                np.dot(stiff_t, [u,u])
+                - sym.FaceMassOperator()(
+                    flux(sym.int_tpair(u))
+                 + flux(sym.bv_tpair(sym.BTAG_ALL, u, bund))
+
+                        )
+                )
+
+
+    bound_op = bind(discr, sym_operator())
+
+    #u = bind(discr, f_gaussian(sym.nodes(dim)))(queue, t=0)
+    u = bind(discr, flow.volume_interpolant()[0:2])(queue, t=0)
+
+    def rhs(t, u):
+        return bound_op(queue, t=t, u=u)
+
+    final_time = 1
+    dt = 0.001
+    print(dt)
+
+    from grudge.shortcuts import set_up_rk4
+    dt_stepper = set_up_rk4("u", dt, u, rhs)
+
+    from grudge.shortcuts import make_visualizer
+    vis = make_visualizer(discr, vis_order=2*order)
+
+    step = 0
+    vis.write_vtk_file("fld-%04d.vtu" % step,
+            #[("u", event.state_component[0])])
+            [("u", u[0]), ("u2", u[1])]
+                )
+
+    for event in dt_stepper.run(t_end=final_time):
+        if isinstance(event, dt_stepper.StateComputed):
+
+            step += 1
+            print(step)
+
+            vis.write_vtk_file("fld-%04d.vtu" % step,
+                    #[("u", event.state_component[0])])
+                    [("u", event.state_component[0]), ("u2", event.state_component[1])]
+                        )
+
+
+if __name__ == "__main__":
+    main()
diff --git a/examples/advection/var-velocity.py b/examples/advection/var-velocity.py
index 00e4f0c8496dced20a21e808e68986ecec1d5be2..7cf7248966a1d8c61ab7d9d73c8263ca52768b54 100644
--- a/examples/advection/var-velocity.py
+++ b/examples/advection/var-velocity.py
@@ -56,9 +56,9 @@ def main(write_output=True, order=4):
 
     sym_x = sym.nodes(2)
 
-    advec_v = join_fields(-1*sym_x[1], sym_x[0]) / 2
+    advec_v = sym_x/sym_x
 
-    flux_type = "upwind"
+    flux_type = "central"
 
     source_center = np.array([0.1, 0.1])
     source_width = 0.05
@@ -109,6 +109,8 @@ def main(write_output=True, order=4):
     dt = dt_factor * h/order**2
     nsteps = (final_time // dt) + 1
     dt = final_time/nsteps + 1e-15
+    dt = 0.001
+    print(dt)
 
     from grudge.shortcuts import set_up_rk4
     dt_stepper = set_up_rk4("u", dt, u, rhs)
@@ -122,7 +124,7 @@ def main(write_output=True, order=4):
         if isinstance(event, dt_stepper.StateComputed):
 
             step += 1
-            if step % 30 == 0:
+            if step % 2 == 0:
                 print(step)
 
                 vis.write_vtk_file("fld-%04d.vtu" % step,
diff --git a/examples/gas_dynamics/euler/sine-wave.py b/examples/gas_dynamics/euler/sine-wave.py
index 451befdbe9613f733ace9324a8770ce8de0d295f..13aeb5870086280fb10877248abec3ee70f65a77 100644
--- a/examples/gas_dynamics/euler/sine-wave.py
+++ b/examples/gas_dynamics/euler/sine-wave.py
@@ -13,17 +13,10 @@
 # You should have received a copy of the GNU General Public License
 # along with this program.  If not, see <http://www.gnu.org/licenses/>.
 
+import pyopencl as cl  # noqa
+import numpy as np
 
-
-
-from __future__ import division
-from __future__ import absolute_import
-from __future__ import print_function
-import numpy
-import numpy.linalg as la
-
-
-
+from grudge import sym, bind, DGDiscretizationWithBoundaries
 
 class SineWave:
     def __init__(self):
@@ -32,161 +25,115 @@ class SineWave:
         self.prandtl = 0.72
         self.spec_gas_const = 287.1
 
-    def __call__(self, t, x_vec):
-        rho = 2 + numpy.sin(x_vec[0] + x_vec[1] + x_vec[2] - 2 * t)
-        velocity = numpy.array([1, 1, 0])
+    def sym_operator(self):
+        t = sym.var("t", sym.DD_SCALAR)
+        sym_nds = sym.nodes(3)
+        rho = 2 + sym.sin(sym_nds[0] + sym_nds[1] + sym_nds[2] - 2 * t)
+        velocity = np.array([1, 1, 0])
         p = 1
-        e = p/(self.gamma-1) + rho/2 * numpy.dot(velocity, velocity)
+        e = p/(self.gamma-1) + rho/2 * np.dot(velocity, velocity)
         rho_u = rho * velocity[0]
         rho_v = rho * velocity[1]
         rho_w = rho * velocity[2]
 
-        from grudge.tools import join_fields
+        from pytools.obj_array import join_fields
         return join_fields(rho, e, rho_u, rho_v, rho_w)
 
     def properties(self):
         return(self.gamma, self.mu, self.prandtl, self.spec_gas_const)
 
-    def volume_interpolant(self, t, discr):
-        return discr.convert_volume(
-                        self(t, discr.nodes.T),
-                        kind=discr.compute_kind)
-
-    def boundary_interpolant(self, t, discr, tag):
-        return discr.convert_boundary(
-                        self(t, discr.get_boundary(tag).nodes.T),
-                         tag=tag, kind=discr.compute_kind)
-
-
-
+def main(order = 3, visualize=False):
+    cl_ctx = cl.create_some_context()
+    queue = cl.CommandQueue(cl_ctx)
 
-def main(final_time=1, write_output=False):
-    from grudge.backends import guess_run_context
-    rcon = guess_run_context()
 
-    from grudge.tools import EOCRecorder, to_obj_array
+    from pytools.convergence import EOCRecorder
     eoc_rec = EOCRecorder()
 
-    if rcon.is_head_rank:
-        from grudge.mesh import make_box_mesh
-        mesh = make_box_mesh((0,0,0), (10,10,10), max_volume=0.5)
-        mesh_data = rcon.distribute_mesh(mesh)
-    else:
-        mesh_data = rcon.receive_mesh()
-
-    for order in [3, 4, 5]:
-        discr = rcon.make_discretization(mesh_data, order=order,
-                        default_scalar_type=numpy.float64)
-
-        from grudge.visualization import SiloVisualizer, VtkVisualizer
-        vis = VtkVisualizer(discr, rcon, "sinewave-%d" % order)
-        #vis = SiloVisualizer(discr, rcon)
-
-        sinewave = SineWave()
-        fields = sinewave.volume_interpolant(0, discr)
-        gamma, mu, prandtl, spec_gas_const = sinewave.properties()
-
-        from grudge.mesh import BTAG_ALL
-        from grudge.models.gas_dynamics import GasDynamicsOperator
-        op = GasDynamicsOperator(dimensions=mesh.dimensions, gamma=gamma, mu=mu,
-                prandtl=prandtl, spec_gas_const=spec_gas_const,
-                bc_inflow=sinewave, bc_outflow=sinewave, bc_noslip=sinewave,
-                inflow_tag=BTAG_ALL, source=None)
-
-        euler_ex = op.bind(discr)
-
-        max_eigval = [0]
-        def rhs(t, q):
-            ode_rhs, speed = euler_ex(t, q)
-            max_eigval[0] = speed
-            return ode_rhs
-        rhs(0, fields)
-
-        if rcon.is_head_rank:
-            print("---------------------------------------------")
-            print("order %d" % order)
-            print("---------------------------------------------")
-            print("#elements=", len(mesh.elements))
-
-        from grudge.timestep import RK4TimeStepper
-        stepper = RK4TimeStepper()
-
-        # diagnostics setup ---------------------------------------------------
-        from pytools.log import LogManager, add_general_quantities, \
-                add_simulation_quantities, add_run_info
-
-        if write_output:
-            log_name = ("euler-sinewave-%(order)d-%(els)d.dat"
-                    % {"order":order, "els":len(mesh.elements)})
-        else:
-            log_name = False
-        logmgr = LogManager(log_name, "w", rcon.communicator)
-        add_run_info(logmgr)
-        add_general_quantities(logmgr)
-        add_simulation_quantities(logmgr)
-        discr.add_instrumentation(logmgr)
-        stepper.add_instrumentation(logmgr)
-
-        logmgr.add_watches(["step.max", "t_sim.max", "t_step.max"])
-
-        # timestep loop -------------------------------------------------------
-        try:
-            from grudge.timestep import times_and_steps
-            step_it = times_and_steps(
-                    final_time=final_time, logmgr=logmgr,
-                    max_dt_getter=lambda t: op.estimate_timestep(discr,
-                        stepper=stepper, t=t, max_eigenvalue=max_eigval[0]))
-
-            for step, t, dt in step_it:
-                #if step % 10 == 0:
-                if write_output:
-                    visf = vis.make_file("sinewave-%d-%04d" % (order, step))
-
-                    #from pyvisfile.silo import DB_VARTYPE_VECTOR
-                    vis.add_data(visf,
-                            [
-                                ("rho", discr.convert_volume(op.rho(fields), kind="numpy")),
-                                ("e", discr.convert_volume(op.e(fields), kind="numpy")),
-                                ("rho_u", discr.convert_volume(op.rho_u(fields), kind="numpy")),
-                                ("u", discr.convert_volume(op.u(fields), kind="numpy")),
-
-                                #("true_rho", op.rho(true_fields)),
-                                #("true_e", op.e(true_fields)),
-                                #("true_rho_u", op.rho_u(true_fields)),
-                                #("true_u", op.u(true_fields)),
-
-                                #("rhs_rho", op.rho(rhs_fields)),
-                                #("rhs_e", op.e(rhs_fields)),
-                                #("rhs_rho_u", op.rho_u(rhs_fields)),
-                                ],
-                            #expressions=[
-                                #("diff_rho", "rho-true_rho"),
-                                #("diff_e", "e-true_e"),
-                                #("diff_rho_u", "rho_u-true_rho_u", DB_VARTYPE_VECTOR),
-
-                                #("p", "0.4*(e- 0.5*(rho_u*u))"),
-                                #],
-                            time=t, step=step
-                            )
-                    visf.close()
-
-                fields = stepper(fields, t, dt, rhs)
-
-        finally:
-            vis.close()
-            logmgr.close()
-            discr.close()
-
-        true_fields = sinewave.volume_interpolant(t, discr)
-        eoc_rec.add_data_point(order, discr.norm(fields-true_fields))
-        print()
-        print(eoc_rec.pretty_print("P.Deg.", "L2 Error"))
+    n = 3
+    from meshmode.mesh.generation import generate_regular_rect_mesh
+    mesh = generate_regular_rect_mesh(a=(-0.5, -0.5, -0.5), b=(0.5, 0.5, 0.5),
+        n=(n, n, n), order=order)
+
+    from meshmode.discretization.poly_element import QuadratureSimplexGroupFactory  # noqa
+    discr = DGDiscretizationWithBoundaries(cl_ctx, mesh, order=order,
+            quad_tag_to_group_factory={
+                #"product": None,
+                "product": QuadratureSimplexGroupFactory(order=3*order)
+                })
+
+
+    sinewave = SineWave()
+    bound_sine = bind(discr, sinewave.sym_operator())
+    fields = bound_sine(queue, t=0)
+    gamma, mu, prandtl, spec_gas_const = sinewave.properties()
+
+    from meshmode.mesh import BTAG_ALL
+    from grudge.models.gas_dynamics import GasDynamicsOperator
+    op = GasDynamicsOperator(dimensions=3, gamma=gamma, mu=mu,
+            prandtl=prandtl, spec_gas_const=spec_gas_const,
+            bc_inflow=sinewave, bc_outflow=sinewave, bc_noslip=sinewave,
+            inflow_tag=BTAG_ALL, source=None)
+    print(sym.pretty(op.sym_operator()))
+
+    euler_ex = bind(discr, op.sym_operator())
+    print(euler_ex)
+
+    max_eigval = [0]
+    def rhs(t, q):
+        #ode_rhs, speed = euler_ex(queue, t=t, q=q)
+        things = euler_ex(queue, t=t, q=q)
+        max_eigval[0] = things[-1]
+        return things[0:-1]
+    rhs(0, fields)
+
+    from grudge.shortcuts import set_up_rk4
+    dt = 0.000001
+    final_t = 0.0001
+    dt_stepper = set_up_rk4("q", dt, fields, rhs)
+
+    if visualize:
+        from grudge.shortcuts import make_visualizer
+        vis = make_visualizer(discr, vis_order=order)
+
+    step = 0
+    for event in dt_stepper.run(t_end=final_t):
+        if isinstance(event, dt_stepper.StateComputed):
+            #assert event.component_id == "u"
+            print(event.component_id)
+            fields = event.state_component
+
+            if visualize:
+                vis.write_vtk_file("fld-sinewave-%d-%04d.vtu" % (order, step), 
+                    [
+                        ("rho", op.rho(fields)),
+                        ("e", op.e(fields)),
+                        ("rho_u", op.rho_u(fields)),
+                        ("u", op.u(fields)),
+
+                        #("true_rho", op.rho(true_fields)),
+                        #("true_e", op.e(true_fields)),
+                        #("true_rho_u", op.rho_u(true_fields)),
+                        #("true_u", op.u(true_fields)),
+
+                        #("rhs_rho", op.rho(rhs_fields)),
+                        #("rhs_e", op.e(rhs_fields)),
+                        #("rhs_rho_u", op.rho_u(rhs_fields)),
+                    ])
+
+            step += 1
+            print(step)
+
+    true_fields = bound_sine(queue, t=t)
+    eoc_rec.add_data_point(order, discr.norm(fields-true_fields))
+    print()
+    print(eoc_rec.pretty_print("P.Deg.", "L2 Error"))
 
 
 
 
 if __name__ == "__main__":
-    main()
+    main(order=2,visualize=True)
 
 
 
diff --git a/examples/gas_dynamics/euler/test.py b/examples/gas_dynamics/euler/test.py
new file mode 100644
index 0000000000000000000000000000000000000000..44d478f12a87481a6a6d5790f47c5a9a620c6889
--- /dev/null
+++ b/examples/gas_dynamics/euler/test.py
@@ -0,0 +1,106 @@
+import numpy as np
+import pyopencl as cl  # noqa
+
+from grudge import sym, bind, DGDiscretizationWithBoundaries
+
+import numpy.linalg as la
+
+
+def test_convergence_maxwell(order, ns, visualize=False):
+    """Test whether 3D maxwells actually converges"""
+
+    cl_ctx = cl.create_some_context()
+    queue = cl.CommandQueue(cl_ctx)
+
+    from pytools.convergence import EOCRecorder
+    eoc_rec = EOCRecorder()
+
+    dims = 2
+    for n in ns:
+        h = 1.0/n
+        dt_factor = 4
+
+        #from meshmode.mesh.generation import generate_regular_rect_mesh
+        #mesh = generate_regular_rect_mesh(a=(-0.5, -0.5), b=(0.5, 0.5),
+            #n=(n, n), order=order)
+
+        from meshmode.mesh.generation import generate_warped_rect_mesh
+        mesh = generate_warped_rect_mesh(dims, order, n)
+
+        discr = DGDiscretizationWithBoundaries(cl_ctx, mesh, order=order)
+
+        advec_v = np.array([0.1, 0.1])
+        norm_v = la.norm(advec_v)
+
+        def f(x):
+            return sym.sin(3*x)
+
+        def u_analytic(x):
+            return f(-np.dot(advec_v, x)/norm_v+sym.var("t", sym.DD_SCALAR)*norm_v)
+
+        from grudge.models.advection import WeakAdvectionOperator
+        analytic_sol = bind(discr, u_analytic(sym.nodes(dims)))
+
+        jac_tag = sym.area_element(2, 2, dd=sym.DOFDesc("vol"))
+        analytic_init = bind(discr, u_analytic(sym.nodes(dims)) / sym.sqrt(jac_tag))
+        J_mul = bind(discr, sym.var("u", sym.DOFDesc("vol")) * sym.sqrt(jac_tag))
+
+        fields = analytic_init(queue, t=0)
+
+        jac_bound = sym.interp(sym.DOFDesc("vol"), sym.BTAG_ALL)(jac_tag)
+        op = WeakAdvectionOperator(advec_v,
+            inflow_u=u_analytic(sym.nodes(dims, sym.BTAG_ALL) / sym.sqrt(jac_bound)),
+            flux_type="central")
+
+        bound_op = bind(discr, op.sym_operator())
+
+        def rhs(t, u):
+            return bound_op(queue, t=t, u=u)
+
+        final_t = 0.1
+        if visualize:
+            final_t = final_t * 10
+
+        dt = dt_factor * h/order**2
+        nsteps = (final_t // dt) + 1
+        dt = final_t/nsteps + 1e-15
+
+        from grudge.shortcuts import set_up_rk4
+        dt_stepper = set_up_rk4("u", dt, fields, rhs)
+
+        print("dt=%g nsteps=%d" % (dt, nsteps))
+
+        norm = bind(discr, sym.norm(2, sym.var("u")))
+
+        step = 0
+        for event in dt_stepper.run(t_end=final_t):
+            if isinstance(event, dt_stepper.StateComputed):
+                assert event.component_id == "u"
+                esc = event.state_component
+
+                step += 1
+                print(step)
+        esc = J_mul(queue, u=esc)
+        
+
+        sol = analytic_sol(queue, t=step * dt)
+        total_error = norm(queue, u=(esc - sol)) / norm(queue, u=sol)
+        eoc_rec.add_data_point(1.0/n, total_error)
+
+        if visualize:
+            from grudge.shortcuts import make_visualizer
+            vis = make_visualizer(discr, vis_order=order)
+            vis.write_vtk_file("fld-thing.vtu", [  ("a", esc), ("s", sol), ("e", (esc-sol))  ])
+            return
+
+
+
+    print(eoc_rec.pretty_print(abscissa_label="h",
+            error_label="L2 Error"))
+
+    #assert eoc_rec.order_estimate() > order
+
+
+if __name__ == "__main__":
+    test_convergence_maxwell(3, [15,20], True)
+
diff --git a/examples/gas_dynamics/euler/vortex.py b/examples/gas_dynamics/euler/vortex.py
index 9f9703f4d77f4c13dc07db0a52eb1ee36fd0ea5d..dd11a51639e76fd92f29931a62ed67a321c36ee7 100644
--- a/examples/gas_dynamics/euler/vortex.py
+++ b/examples/gas_dynamics/euler/vortex.py
@@ -21,181 +21,141 @@ from __future__ import absolute_import
 from __future__ import print_function
 import numpy
 
+import pyopencl as cl  # noqa
 
+from grudge import sym, bind, DGDiscretizationWithBoundaries
 
-
-def main(write_output=True):
+def main(write_output=True, visualize=False):
     from pytools import add_python_path_relative_to_script
     add_python_path_relative_to_script("..")
+    cl_ctx = cl.create_some_context()
+    queue = cl.CommandQueue(cl_ctx)
 
-    from grudge.backends import guess_run_context
-    rcon = guess_run_context()
 
-    from grudge.tools import EOCRecorder
+    from pytools.convergence import EOCRecorder
     eoc_rec = EOCRecorder()
 
-    if rcon.is_head_rank:
-        from grudge.mesh.generator import \
-                make_rect_mesh, \
-                make_centered_regular_rect_mesh
+    from meshmode.mesh.generation import generate_regular_rect_mesh
 
-        refine = 4
-        mesh = make_centered_regular_rect_mesh((0,-5), (10,5), n=(9,9),
-                post_refine_factor=refine)
-        mesh_data = rcon.distribute_mesh(mesh)
-    else:
-        mesh_data = rcon.receive_mesh()
 
     for order in [3, 4, 5]:
+        refine = 4
+        mesh = generate_regular_rect_mesh(a=(0,-5), b=(10,5), n=(6,6), order=order)
         from gas_dynamics_initials import Vortex
         flow = Vortex()
 
         from grudge.models.gas_dynamics import (
                 GasDynamicsOperator, PolytropeEOS, GammaLawEOS)
 
-        from grudge.mesh import BTAG_ALL
+        from meshmode.mesh import BTAG_ALL
+        from meshmode.discretization.poly_element import QuadratureSimplexGroupFactory  
         # works equally well for GammaLawEOS
         op = GasDynamicsOperator(dimensions=2, mu=flow.mu,
                 prandtl=flow.prandtl, spec_gas_const=flow.spec_gas_const,
                 equation_of_state=PolytropeEOS(flow.gamma),
-                bc_inflow=flow, bc_outflow=flow, bc_noslip=flow,
+                bc_inflow=flow.volume_interpolant(), bc_outflow=flow.volume_interpolant(), bc_noslip=flow.volume_interpolant(),
                 inflow_tag=BTAG_ALL, source=None)
 
-        discr = rcon.make_discretization(mesh_data, order=order,
-                        default_scalar_type=numpy.float64,
-                        quad_min_degrees={
-                            "gasdyn_vol": 3*order,
-                            "gasdyn_face": 3*order,
-                            },
-                        tune_for=op.sym_operator(),
-                        debug=["cuda_no_plan"])
+        discr = DGDiscretizationWithBoundaries(cl_ctx,mesh, order=order,
+                        quad_tag_to_group_factory={
+                "product": QuadratureSimplexGroupFactory(order=3*order)
+                })
+
+        if visualize:
+            from grudge.shortcuts import make_visualizer
+            vis = make_visualizer(discr, vis_order=order)
 
-        from grudge.visualization import SiloVisualizer, VtkVisualizer
-        vis = VtkVisualizer(discr, rcon, "vortex-%d" % order)
-        #vis = SiloVisualizer(discr, rcon)
 
-        fields = flow.volume_interpolant(0, discr)
+        bound_flow_vol = bind(discr, flow.volume_interpolant())
+        fields = bound_flow_vol(queue, t=1)
+        print(fields)
 
-        euler_ex = op.bind(discr)
+        euler_ex = bind(discr, op.sym_operator())
 
         max_eigval = [0]
         def rhs(t, q):
-            ode_rhs, speed = euler_ex(t, q)
+            #ode_rhs, speed = euler_ex(queue, t=t, q=q)
+            ting = euler_ex(queue, t=t, q=q)
+            ode_rhs = ting[:-1]
+            speed = ting[-1]
+            #pu.db
             max_eigval[0] = speed
             return ode_rhs
         rhs(0, fields)
 
-        if rcon.is_head_rank:
-            print("---------------------------------------------")
-            print("order %d" % order)
-            print("---------------------------------------------")
-            print("#elements=", len(mesh.elements))
-
 
         # limiter ------------------------------------------------------------
-        from grudge.models.gas_dynamics import SlopeLimiter1NEuler
-        limiter = SlopeLimiter1NEuler(discr, flow.gamma, 2, op)
+        #from grudge.models.gas_dynamics import SlopeLimiter1NEuler
+        #limiter = SlopeLimiter1NEuler(discr, flow.gamma, 2, op)
 
-        from grudge.timestep.runge_kutta import SSP3TimeStepper
+        #from grudge.timestep.runge_kutta import SSP3TimeStepper
         #stepper = SSP3TimeStepper(limiter=limiter)
-        stepper = SSP3TimeStepper(
-                vector_primitive_factory=discr.get_vector_primitive_factory())
-
-        #from grudge.timestep import RK4TimeStepper
-        #stepper = RK4TimeStepper()
-
-        # diagnostics setup ---------------------------------------------------
-        from pytools.log import LogManager, add_general_quantities, \
-                add_simulation_quantities, add_run_info
-
-        if write_output:
-            log_file_name = "euler-%d.dat" % order
-        else:
-            log_file_name = None
-
-        logmgr = LogManager(log_file_name, "w", rcon.communicator)
-        add_run_info(logmgr)
-        add_general_quantities(logmgr)
-        add_simulation_quantities(logmgr)
-        discr.add_instrumentation(logmgr)
-        stepper.add_instrumentation(logmgr)
-
-        logmgr.add_watches(["step.max", "t_sim.max", "t_step.max"])
-
-        # timestep loop -------------------------------------------------------
-        try:
-            final_time = flow.final_time
-            from grudge.timestep import times_and_steps
-            step_it = times_and_steps(
-                    final_time=final_time, logmgr=logmgr,
-                    max_dt_getter=lambda t: op.estimate_timestep(discr,
-                        stepper=stepper, t=t, max_eigenvalue=max_eigval[0]))
-
-            print("run until t=%g" % final_time)
-            for step, t, dt in step_it:
-                if step % 10 == 0 and write_output:
-                #if False:
-                    visf = vis.make_file("vortex-%d-%04d" % (order, step))
-
-                    #true_fields = vortex.volume_interpolant(t, discr)
-
-                    from pyvisfile.silo import DB_VARTYPE_VECTOR
-                    vis.add_data(visf,
-                            [
-                                ("rho", discr.convert_volume(op.rho(fields), kind="numpy")),
-                                ("e", discr.convert_volume(op.e(fields), kind="numpy")),
-                                ("rho_u", discr.convert_volume(op.rho_u(fields), kind="numpy")),
-                                ("u", discr.convert_volume(op.u(fields), kind="numpy")),
-
-                                #("true_rho", discr.convert_volume(op.rho(true_fields), kind="numpy")),
-                                #("true_e", discr.convert_volume(op.e(true_fields), kind="numpy")),
-                                #("true_rho_u", discr.convert_volume(op.rho_u(true_fields), kind="numpy")),
-                                #("true_u", discr.convert_volume(op.u(true_fields), kind="numpy")),
-
-                                #("rhs_rho", discr.convert_volume(op.rho(rhs_fields), kind="numpy")),
-                                #("rhs_e", discr.convert_volume(op.e(rhs_fields), kind="numpy")),
-                                #("rhs_rho_u", discr.convert_volume(op.rho_u(rhs_fields), kind="numpy")),
-                                ],
-                            #expressions=[
-                                #("diff_rho", "rho-true_rho"),
-                                #("diff_e", "e-true_e"),
-                                #("diff_rho_u", "rho_u-true_rho_u", DB_VARTYPE_VECTOR),
-
-                                #("p", "0.4*(e- 0.5*(rho_u*u))"),
-                                #],
-                            time=t, step=step
-                            )
-                    visf.close()
-
-                fields = stepper(fields, t, dt, rhs)
-                #fields = limiter(fields)
-
-                assert not numpy.isnan(numpy.sum(fields[0]))
-
-            true_fields = flow.volume_interpolant(final_time, discr)
-            l2_error = discr.norm(fields-true_fields)
-            l2_error_rho = discr.norm(op.rho(fields)-op.rho(true_fields))
-            l2_error_e = discr.norm(op.e(fields)-op.e(true_fields))
-            l2_error_rhou = discr.norm(op.rho_u(fields)-op.rho_u(true_fields))
-            l2_error_u = discr.norm(op.u(fields)-op.u(true_fields))
-
-            eoc_rec.add_data_point(order, l2_error)
-            print()
-            print(eoc_rec.pretty_print("P.Deg.", "L2 Error"))
-
-            logmgr.set_constant("l2_error", l2_error)
-            logmgr.set_constant("l2_error_rho", l2_error_rho)
-            logmgr.set_constant("l2_error_e", l2_error_e)
-            logmgr.set_constant("l2_error_rhou", l2_error_rhou)
-            logmgr.set_constant("l2_error_u", l2_error_u)
-            logmgr.set_constant("refinement", refine)
-
-        finally:
-            if write_output:
-                vis.close()
-
-            logmgr.close()
-            discr.close()
+        #stepper = SSP3TimeStepper(
+                #vector_primitive_factory=discr.get_vector_primitive_factory())
+        from grudge.shortcuts import set_up_rk4
+        dt = 0.001
+        final_t = 0.1
+        dt_stepper = set_up_rk4("q", dt, fields, rhs)
+
+        if visualize:
+            from grudge.shortcuts import make_visualizer
+            vis = make_visualizer(discr, vis_order=order)
+
+        norm = bind(discr, sym.norm(2, sym.var("q")))
+
+        vis.write_vtk_file("fld-sinewave-%d-%04d.vtu" % (order, 0),
+            [
+                ("rho", op.rho(fields)),
+                ("e", op.e(fields)),
+                ("rho_u", op.rho_u(fields)),
+                ("u", op.u(fields))
+                #("true_rho", op.rho(true_fields)),
+                #("true_e", op.e(true_fields)),
+                #("true_rho_u", op.rho_u(true_fields)),
+                #("true_u", op.u(true_fields)
+                #("rhs_rho", op.rho(rhs_fields)),
+                #("rhs_e", op.e(rhs_fields)),
+                #("rhs_rho_u", op.rho_u(rhs_fields)),
+            ])
+        step = 1
+        for event in dt_stepper.run(t_end=final_t):
+            if isinstance(event, dt_stepper.StateComputed):
+                assert event.component_id == "q"
+                fields = event.state_component
+
+                if visualize:
+                    vis.write_vtk_file("fld-sinewave-%d-%04d.vtu" % (order, step),
+                        [
+                            ("rho", op.rho(fields)),
+                            ("e", op.e(fields)),
+                            ("rho_u", op.rho_u(fields)),
+                            ("u", op.u(fields)),
+
+                            #("true_rho", op.rho(true_fields)),
+                            #("true_e", op.e(true_fields)),
+                            #("true_rho_u", op.rho_u(true_fields)),
+                            #("true_u", op.u(true_fields)),
+
+                            #("rhs_rho", op.rho(rhs_fields)),
+                            #("rhs_e", op.e(rhs_fields)),
+                            #("rhs_rho_u", op.rho_u(rhs_fields)),
+                        ])
+                step += 1
+                print(step)
+
+
+
+                #true_fields = bound_flow_vol(queue, t=final_t)
+                #l2_error = norm(queue, q = fields-true_fields)
+                #l2_error_rho = discr.norm(op.rho(fields)-op.rho(true_fields))
+                #l2_error_e = discr.norm(op.e(fields)-op.e(true_fields))
+                #l2_error_rhou = discr.norm(op.rho_u(fields)-op.rho_u(true_fields))
+                #l2_error_u = discr.norm(op.u(fields)-op.u(true_fields))
+
+                #eoc_rec.add_data_point(order, l2_error)
+                #print(eoc_rec.pretty_print("P.Deg.", "L2 Error"))
+
+
 
     # after order loop
     assert eoc_rec.estimate_order_of_convergence()[0,1] > 6
@@ -204,7 +164,7 @@ def main(write_output=True):
 
 
 if __name__ == "__main__":
-    main()
+    main(visualize=True)
 
 
 
diff --git a/examples/gas_dynamics/euler/vortex2.py b/examples/gas_dynamics/euler/vortex2.py
new file mode 100644
index 0000000000000000000000000000000000000000..81f19c296cea3ef93a2ec50ce5f2a41aa4a04289
--- /dev/null
+++ b/examples/gas_dynamics/euler/vortex2.py
@@ -0,0 +1,154 @@
+# Copyright (C) 2008 Andreas Kloeckner
+#
+# This program is free software: you can redistribute it and/or modify
+# it under the terms of the GNU General Public License as published by
+# the Free Software Foundation, either version 3 of the License, or
+# (at your option) any later version.
+#
+# This program is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+# GNU General Public License for more details.
+#
+# You should have received a copy of the GNU General Public License
+# along with this program.  If not, see <http://www.gnu.org/licenses/>.
+
+
+
+from __future__ import division
+from __future__ import absolute_import
+from __future__ import print_function
+import numpy
+
+import pyopencl as cl  # noqa
+
+from grudge import sym, bind, DGDiscretizationWithBoundaries
+
+def main(write_output=True, visualize=False):
+    from pytools import add_python_path_relative_to_script
+    add_python_path_relative_to_script("..")
+    cl_ctx = cl.create_some_context()
+    queue = cl.CommandQueue(cl_ctx)
+
+
+    from pytools.convergence import EOCRecorder
+    eoc_rec = EOCRecorder()
+
+    from meshmode.mesh.generation import generate_regular_rect_mesh
+
+
+    for order in [5]:
+        mesh = generate_regular_rect_mesh(a=(0,-5), b=(10,5), n=(5,5), order=order)
+        from gas_dynamics_initials import Vortex
+        flow = Vortex()  # NOT ACTUALLY A VORTEX, constant velocity
+
+        from grudge.models.euler2 import EulerOperator
+        from grudge.models.gas_dynamics import (
+                GasDynamicsOperator, PolytropeEOS, GammaLawEOS)
+
+        from meshmode.mesh import BTAG_ALL
+        from meshmode.discretization.poly_element import QuadratureSimplexGroupFactory  
+        # works equally well for GammaLawEOS
+        op = EulerOperator(dimensions=2, mu=flow.mu,
+                
+                equation_of_state=PolytropeEOS(flow.gamma),
+                bc_inflow=flow.volume_interpolant(),
+                inflow_tag=BTAG_ALL)
+
+        discr = DGDiscretizationWithBoundaries(cl_ctx,mesh, order=order,
+                        quad_tag_to_group_factory={
+                "product": QuadratureSimplexGroupFactory(order=3*order)
+                })
+
+        if visualize:
+            from grudge.shortcuts import make_visualizer
+            vis = make_visualizer(discr, vis_order=order)
+
+
+        bound_flow_vol = bind(discr, flow.volume_interpolant())
+        fields = bound_flow_vol(queue, t=0)
+        print(fields)
+
+        euler_ex = bind(discr, op.sym_operator())
+
+        max_eigval = [0]
+        def rhs(t, q):
+            #ode_rhs, speed = euler_ex(queue, t=t, q=q)
+            return euler_ex(queue, t=t, q=q)
+            ode_rhs = ting[:-1]
+            speed = ting[-1]
+            #pu.db
+            max_eigval[0] = speed
+            return ode_rhs
+        rhs(0, fields)
+
+        from grudge.shortcuts import set_up_rk4
+        dt = 0.01
+        final_t = 2
+        dt_stepper = set_up_rk4("q", dt, fields, rhs)
+
+        if visualize:
+            from grudge.shortcuts import make_visualizer
+            vis = make_visualizer(discr, vis_order=order)
+
+        norm = bind(discr, sym.norm(2, sym.var("q")))
+
+        vis.write_vtk_file("fld-sinewave-%d-%04d.vtu" % (order, 0),
+            [
+                ("rho", op.rho(fields)),
+                ("e", op.e(fields)),
+                ("rho_u", op.rho_u(fields)),
+                ("u", op.u(fields))
+                #("true_rho", op.rho(true_fields)),
+                #("true_e", op.e(true_fields)),
+                #("true_rho_u", op.rho_u(true_fields)),
+                #("true_u", op.u(true_fields)
+                #("rhs_rho", op.rho(rhs_fields)),
+                #("rhs_e", op.e(rhs_fields)),
+                #("rhs_rho_u", op.rho_u(rhs_fields)),
+            ])
+        step = 1
+        for event in dt_stepper.run(t_end=final_t):
+            if isinstance(event, dt_stepper.StateComputed):
+                assert event.component_id == "q"
+                fields = event.state_component
+
+                if visualize:
+                    vis.write_vtk_file("fld-sinewave-%d-%04d.vtu" % (order, step),
+                        [
+                            ("rho", op.rho(fields)),
+                            ("e", op.e(fields)),
+                            ("rho_u", op.rho_u(fields)),
+                            ("u", op.u(fields)),
+
+                            #("true_rho", op.rho(true_fields)),
+                            #("true_e", op.e(true_fields)),
+                            #("true_rho_u", op.rho_u(true_fields)),
+                            #("true_u", op.u(true_fields)),
+
+                            #("rhs_rho", op.rho(rhs_fields)),
+                            #("rhs_e", op.e(rhs_fields)),
+                            #("rhs_rho_u", op.rho_u(rhs_fields)),
+                        ])
+                step += 1
+                print(step)
+
+
+
+                #true_fields = bound_flow_vol(queue, t=final_t)
+                #l2_error = norm(queue, q = fields-true_fields)
+                #l2_error_rho = discr.norm(op.rho(fields)-op.rho(true_fields))
+                #l2_error_e = discr.norm(op.e(fields)-op.e(true_fields))
+                #l2_error_rhou = discr.norm(op.rho_u(fields)-op.rho_u(true_fields))
+                #l2_error_u = discr.norm(op.u(fields)-op.u(true_fields))
+
+                #eoc_rec.add_data_point(order, l2_error)
+                #print(eoc_rec.pretty_print("P.Deg.", "L2 Error"))
+
+
+
+    # after order loop
+    assert eoc_rec.estimate_order_of_convergence()[0,1] > 6
+
+if __name__ == "__main__":
+    main(visualize=True)
diff --git a/examples/gas_dynamics/gas_dynamics_initials.py b/examples/gas_dynamics/gas_dynamics_initials.py
index 05ab91a314415a80d8a30452e2d68c3e26d6d378..98cf671ba37ef21e672db095e6e04a8b3a7dfaae 100644
--- a/examples/gas_dynamics/gas_dynamics_initials.py
+++ b/examples/gas_dynamics/gas_dynamics_initials.py
@@ -21,6 +21,7 @@ from __future__ import absolute_import
 import numpy
 import numpy.linalg as la
 from six.moves import range
+from grudge import sym
 
 
 
@@ -120,83 +121,39 @@ class UniformMachFlow:
                          tag=tag, kind=discr.compute_kind,
                          dtype=discr.default_scalar_type)
 
-class Vortex:
-    def __init__(self):
-        self.beta = 5
-        self.gamma = 1.4
-        self.center = numpy.array([5, 0])
-        self.velocity = numpy.array([1, 0])
 
-        self.mu = 0
-        self.prandtl = 0.72
-        self.spec_gas_const = 287.1
 
-    def __call__(self, t, x_vec):
-        vortex_loc = self.center + t*self.velocity
-
-        # coordinates relative to vortex center
-        x_rel = x_vec[0] - vortex_loc[0]
-        y_rel = x_vec[1] - vortex_loc[1]
 
-        # Y.C. Zhou, G.W. Wei / Journal of Computational Physics 189 (2003) 159
-        # also JSH/TW Nodal DG Methods, p. 209
 
-        from math import pi
-        r = numpy.sqrt(x_rel**2+y_rel**2)
-        expterm = self.beta*numpy.exp(1-r**2)
-        u = self.velocity[0] - expterm*y_rel/(2*pi)
-        v = self.velocity[1] + expterm*x_rel/(2*pi)
-        rho = (1-(self.gamma-1)/(16*self.gamma*pi**2)*expterm**2)**(1/(self.gamma-1))
-        p = rho**self.gamma
 
-        e = p/(self.gamma-1) + rho/2*(u**2+v**2)
 
-        from grudge.tools import join_fields
-        return join_fields(rho, e, rho*u, rho*v)
-
-    def volume_interpolant(self, t, discr):
-        return discr.convert_volume(
-                        self(t, discr.nodes.T
-                            .astype(discr.default_scalar_type)),
-                        kind=discr.compute_kind)
-
-    def boundary_interpolant(self, t, discr, tag):
-        return discr.convert_boundary(
-                        self(t, discr.get_boundary(tag).nodes.T
-                            .astype(discr.default_scalar_type)),
-                         tag=tag, kind=discr.compute_kind)
-
-
-
-
-
-
-
-class Vortex:
+class Vortex: # NOT ACTUALLY A VORTEX, constant velocity
     def __init__(self):
         self.beta = 5
         self.gamma = 1.4
         self.center = numpy.array([5, 0])
-        self.velocity = numpy.array([1, 0])
+        self.velocity = numpy.array([0.3, 0.3])
         self.final_time = 0.5
 
         self.mu = 0
         self.prandtl = 0.72
         self.spec_gas_const = 287.1
 
-    def __call__(self, t, x_vec):
+    def __call__(self):
+        x_vec = sym.nodes(2)
+        t = sym.var("t", sym.DD_SCALAR)
+        from pytools.obj_array import join_fields
         vortex_loc = self.center + t*self.velocity
 
         # coordinates relative to vortex center
         x_rel = x_vec[0] - vortex_loc[0]
         y_rel = x_vec[1] - vortex_loc[1]
-
-        # Y.C. Zhou, G.W. Wei / Journal of Computational Physics 189 (2003) 159
+# Y.C. Zhou, G.W. Wei / Journal of Computational Physics 189 (2003) 159
         # also JSH/TW Nodal DG Methods, p. 209
 
         from math import pi
-        r = numpy.sqrt(x_rel**2+y_rel**2)
-        expterm = self.beta*numpy.exp(1-r**2)
+        r = sym.sqrt(x_rel**2+y_rel**2)
+        expterm = self.beta*sym.exp(1-r**2)
         u = self.velocity[0] - expterm*y_rel/(2*pi)
         v = self.velocity[1] + expterm*x_rel/(2*pi)
         rho = (1-(self.gamma-1)/(16*self.gamma*pi**2)*expterm**2)**(1/(self.gamma-1))
@@ -204,20 +161,10 @@ class Vortex:
 
         e = p/(self.gamma-1) + rho/2*(u**2+v**2)
 
-        from grudge.tools import join_fields
         return join_fields(rho, e, rho*u, rho*v)
 
-    def volume_interpolant(self, t, discr):
-        return discr.convert_volume(
-                        self(t, discr.nodes.T
-                            .astype(discr.default_scalar_type)),
-                        kind=discr.compute_kind)
-
-    def boundary_interpolant(self, t, discr, tag):
-        return discr.convert_boundary(
-                        self(t, discr.get_boundary(tag).nodes.T
-                            .astype(discr.default_scalar_type)),
-                         tag=tag, kind=discr.compute_kind)
+    def volume_interpolant(self):
+        return self()
 
 
 
diff --git a/examples/heat/heat2.py b/examples/heat/heat2.py
new file mode 100644
index 0000000000000000000000000000000000000000..54306eedae3cb1454e9cb2f8ab9eb679e60acc7e
--- /dev/null
+++ b/examples/heat/heat2.py
@@ -0,0 +1,124 @@
+# Copyright (C) 2007 Andreas Kloeckner
+#
+# This program is free software: you can redistribute it and/or modify
+# it under the terms of the GNU General Public License as published by
+# the Free Software Foundation, either version 3 of the License, or
+# (at your option) any later version.
+#
+# This program is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+# GNU General Public License for more details.
+#
+# You should have received a copy of the GNU General Public License
+# along with this program.  If not, see <http://www.gnu.org/licenses/>.
+
+
+import numpy as np
+import pyopencl as cl  # noqa
+import pyopencl.array  # noqa
+import pyopencl.clmath  # noqa
+
+import pytest  # noqa
+
+from pyopencl.tools import (  # noqa
+                pytest_generate_tests_for_pyopencl as pytest_generate_tests)
+
+import logging
+logger = logging.getLogger(__name__)
+
+from grudge import sym, bind, DGDiscretizationWithBoundaries
+
+import numpy.linalg as la
+
+
+
+
+def main(write_output=True, order=4):
+    cl_ctx = cl.create_some_context()
+    queue = cl.CommandQueue(cl_ctx)
+
+    dim = 2
+
+
+    from meshmode.mesh.generation import generate_regular_rect_mesh
+    mesh = generate_regular_rect_mesh(a=(-0.5, -0.5), b=(0.5, 0.5),
+            n=(20, 20), order=order)
+
+    dt_factor = 4
+    h = 1/20
+
+    discr = DGDiscretizationWithBoundaries(cl_ctx, mesh, order=order)
+
+    c = np.array([0.1,0.1])
+    norm_c = la.norm(c)
+
+
+    flux_type = "central"
+         
+
+    def f(x):
+        return sym.If(
+            sym.Comparison(sym.norm(x), ">", 0.1),
+            1,2)
+
+    def u_analytic(x):
+        return f(x)
+
+    from grudge.models.diffusion import DiffusionOperator
+    from meshmode.mesh import BTAG_ALL, BTAG_NONE
+    
+    discr = DGDiscretizationWithBoundaries(cl_ctx, mesh, order=order)
+    op = DiffusionOperator(2,
+            dirichlet_tag="dirichlet",
+            dirichlet_bc=0,#bad
+            neumann_tag="neumann",
+            neumann_bc=1,#bad
+            )
+
+
+    bound_op = bind(discr, op.m_sym_operator())
+
+    u = bind(discr, u_analytic(sym.nodes(dim)))(queue)
+
+    def rhs(t, u):
+        return bound_op(queue, t=t, u=u)
+
+    final_time = 0.3
+
+    dt = dt_factor * h/order**2
+    nsteps = (final_time // dt) + 1
+    dt = final_time/nsteps + 1e-15
+
+
+    from grudge.shortcuts import set_up_rk4
+    dt_stepper = set_up_rk4("u", dt, u, rhs)
+
+    last_u = None
+
+    from grudge.shortcuts import make_visualizer
+    vis = make_visualizer(discr, vis_order=order)
+
+    step = 0
+
+    norm = bind(discr, sym.norm(2, sym.var("u")))
+
+    for event in dt_stepper.run(t_end=final_time):
+        if isinstance(event, dt_stepper.StateComputed):
+
+            step += 1
+
+            #print(step, event.t, norm(queue, u=event.state_component[0]))
+            vis.write_vtk_file("fld-%04d.vtu" % step,
+                    [  ("u", event.state_component) ])
+
+
+
+
+
+
+
+if __name__ == "__main__":
+    main()
+
+
diff --git a/examples/maxwell/cavities.py b/examples/maxwell/cavities.py
index 85a7c3dec0cef99b4a72fe9c23439c4212c0c32d..a8763e8461c8836c9cf1e04a93af2612e37c0a59 100644
--- a/examples/maxwell/cavities.py
+++ b/examples/maxwell/cavities.py
@@ -28,7 +28,7 @@ THE SOFTWARE.
 import numpy as np
 import pyopencl as cl
 from grudge.shortcuts import set_up_rk4
-from grudge import sym, bind, Discretization
+from grudge import sym, bind, DGDiscretizationWithBoundaries
 
 from grudge.models.em import get_rectangular_cavity_mode
 
@@ -43,7 +43,7 @@ def main(dims, write_output=True, order=4):
             b=(1.0,)*dims,
             n=(5,)*dims)
 
-    discr = Discretization(cl_ctx, mesh, order=order)
+    discr = DGDiscretizationWithBoundaries(cl_ctx, mesh, order=order)
 
     if 0:
         epsilon0 = 8.8541878176e-12  # C**2 / (N m**2)
diff --git a/examples/wave/wave-min.py b/examples/wave/wave-min.py
index 89f4bfa836abe91e305802a455fd19a729cac172..e366a99273b86c1cf25f726dc4fc422700266c8e 100644
--- a/examples/wave/wave-min.py
+++ b/examples/wave/wave-min.py
@@ -35,12 +35,12 @@ def main(write_output=True, order=4):
     cl_ctx = cl.create_some_context()
     queue = cl.CommandQueue(cl_ctx)
 
-    dims = 3
+    dims = 2
     from meshmode.mesh.generation import generate_regular_rect_mesh
     mesh = generate_regular_rect_mesh(
             a=(-0.5,)*dims,
             b=(0.5,)*dims,
-            n=(16,)*dims)
+            n=(9,)*dims)
 
     if mesh.dim == 2:
         dt = 0.04
@@ -114,7 +114,7 @@ def main(write_output=True, order=4):
 
             print(step, event.t, norm(queue, u=event.state_component[0]),
                     time()-t_last_step)
-            if step % 10 == 0:
+            if step % 3 == 0:
                 vis.write_vtk_file("fld-%04d.vtu" % step,
                         [
                             ("u", event.state_component[0]),
@@ -124,4 +124,4 @@ def main(write_output=True, order=4):
 
 
 if __name__ == "__main__":
-    main()
+    main(order = 4)
diff --git a/grudge/execution.py b/grudge/execution.py
index 2da2e508cc8ddf116eb4ffb30ecfd3f8835f957b..5c7228a6e9fddeb1a6710d33b174bb4d26c37acc 100644
--- a/grudge/execution.py
+++ b/grudge/execution.py
@@ -114,6 +114,7 @@ class ExecutionMapper(mappers.Evaluator,
             cached_name += "abs"
         else:
             cached_name += expr.function.name
+        #print(cached_name + ":")
 
         @memoize_in(self.bound_op, cached_name)
         def knl():
@@ -128,6 +129,8 @@ class ExecutionMapper(mappers.Evaluator,
         assert len(args) == 1
         evt, (out,) = knl()(self.queue, a=args[0])
 
+        #print(out)
+
         return out
 
     def map_nodal_sum(self, op, field_expr):
@@ -233,15 +236,72 @@ class ExecutionMapper(mappers.Evaluator,
         return result
 
     def map_elementwise_max(self, op, field_expr):
-        from grudge._internal import perform_elwise_max
         field = self.rec(field_expr)
 
-        out = self.discr.volume_zeros(dtype=field.dtype)
-        for eg in self.discr.element_groups:
-            perform_elwise_max(eg.ranges, field, out)
+        from grudge.tools import is_zero
+        if is_zero(field):
+            return 0
+
+        @memoize_in(self.bound_op, "elwise_max_knl")
+        def knl():
+            knl = lp.make_kernel(
+                """{[k,i,j]:
+                    0<=k<nelements and
+                    0<=i<ndiscr_nodes_out and
+                    0<=j<ndiscr_nodes_in}""",
+                "result[k,i] = max(j, vec[k, j])",
+                default_offset=lp.auto, name="diff")
+
+            knl = lp.split_iname(knl, "i", 16, inner_tag="l.0")
+            return lp.tag_inames(knl, dict(k="g.0"))
+
+        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)
+
+        for in_grp, out_grp in zip(in_discr.groups, out_discr.groups):
+            knl()(self.queue, result=out_grp.view(result),
+                    vec=in_grp.view(field))
+
+        return result
+
+        inp = self.rec(field_expr)
+        print(inp.shape)
+
+        import pymbolic.primitives as p
+        var = p.Variable
+
+        i = var("i")
+        @memoize_in(self.bound_op, "face_elwisemax_knl")
+        def knl():
+            knl = lp.make_kernel(
+                "{[i,j]: 0<=i<n and 0<=j<n}",
+                """
+                    out[i] = max(j,inp)
+                """)
+                #[
+                    #lp.Assignment(var("out")[i],
+                        #lp.Reduction(lp.library.reduction.MaxReductionOperation(), ("j",), var("inp"), True))
+                #])
+            return lp.split_iname(knl, "i", 128, outer_tag="g.0", inner_tag="l.0")  
+
+        evt, (out,) = knl()(self.queue, inp = inp, n=inp.shape[0])
 
         return out
 
+        
+        #from grudge._internal import perform_elwise_max
+        #field = self.rec(field_expr)
+
+        #out = self.discr.volume_zeros(dtype=field.dtype)
+        #for eg in self.discr.element_groups:
+            #perform_elwise_max(eg.ranges, field, out)
+
+        #return out
+
     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)
@@ -330,7 +390,12 @@ class ExecutionMapper(mappers.Evaluator,
                         % (name, type(v)))
 
         kwargs["grdg_n"] = discr.nnodes
+        #kdescr.loopy_kernel = lp.set_options(kdescr.loopy_kernel, trace_assignment_values=True)
         evt, result_dict = kdescr.loopy_kernel(self.queue, **kwargs)
+        #if "expr_217" in result_dict.keys():
+            #knl = lp.set_options(kdescr.loopy_kernel, "write_cl")
+            #knl(self.queue, **kwargs)
+            #pu.db
         return list(result_dict.items()), []
 
     def map_insn_assign(self, insn):
diff --git a/grudge/models/advection.py b/grudge/models/advection.py
index 091b9870c4b82a3fde2fc21b84bf92e5ec62ab17..d496760c189f5035e8012b7a081dbf51474fbae0 100644
--- a/grudge/models/advection.py
+++ b/grudge/models/advection.py
@@ -165,36 +165,34 @@ class VariableCoefficientAdvectionOperator(HyperbolicOperator):
             raise ValueError("invalid flux type")
 
     def sym_operator(self):
+        def mflux(pair):
+            normal = sym.normal(pair.dd, self.ambient_dim)
+            return pair.avg*np.dot(normal, [1,1])
+
+        def flux(pair):
+            return sym.interp(pair.dd, "all_faces")(
+                    mflux(pair))
         u = sym.var("u")
+        stiff_t = sym.stiffness_t(self.ambient_dim)
+        return sym.InverseMassOperator()(
+                np.dot(stiff_t, [u,u])
+                - sym.FaceMassOperator()(
+                    flux(sym.int_tpair(u))
+                ))
 
         # boundary conditions -------------------------------------------------
         bc_in = self.inflow_u
 
-        all_faces_dd = sym.DOFDesc(sym.FACE_RESTR_ALL, self.quad_tag)
         boundary_dd = sym.DOFDesc(sym.BTAG_ALL, self.quad_tag)
 
-        def flux(pair):
-            return sym.interp(pair.dd, all_faces_dd)(
-                    self.flux(pair))
 
         quad_dd = sym.DOFDesc("vol", self.quad_tag)
 
         to_quad = sym.interp("vol", quad_dd)
 
-        stiff_t = sym.stiffness_t(self.ambient_dim, quad_dd, "vol")
         quad_v = to_quad(self.v)
         quad_u = to_quad(u)
 
-        return sym.InverseMassOperator()(
-                (stiff_t[0](quad_u * quad_v[0]) + stiff_t[1](quad_u * quad_v[1]))
-                - sym.FaceMassOperator(all_faces_dd, "vol")(
-                    flux(sym.int_tpair(u, self.quad_tag))
-                    + flux(sym.bv_tpair(boundary_dd, u, bc_in))
-
-                    # FIXME: Add back support for inflow/outflow tags
-                    #+ flux(sym.bv_tpair(self.inflow_tag, u, bc_in))
-                    #+ flux(sym.bv_tpair(self.outflow_tag, u, bc_out))
-                ))
 # }}}
 
 # vim: foldmethod=marker
diff --git a/grudge/models/diffusion.py b/grudge/models/diffusion.py
index 04c5cf57c5368db4b88dfd4008b8962490f4d317..81319a92413e8b0bd36175c71c06f19ea5a731cf 100644
--- a/grudge/models/diffusion.py
+++ b/grudge/models/diffusion.py
@@ -30,18 +30,17 @@ THE SOFTWARE.
 
 import numpy
 
-import grudge.data
 from grudge.models import TimeDependentOperator
 from grudge.models.poisson import LaplacianOperatorBase
-from grudge.second_order import CentralSecondDerivative
+from grudge.models.second_order import CentralSecondDerivative
 
 
 
 
 class DiffusionOperator(TimeDependentOperator, LaplacianOperatorBase):
     def __init__(self, dimensions, diffusion_tensor=None,
-            dirichlet_bc=grudge.data.make_tdep_constant(0), dirichlet_tag="dirichlet",
-            neumann_bc=grudge.data.make_tdep_constant(0), neumann_tag="neumann",
+            dirichlet_bc=0, dirichlet_tag="dirichlet",
+            neumann_bc=0, neumann_tag="neumann",
             scheme=CentralSecondDerivative()):
         self.dimensions = dimensions
 
@@ -61,11 +60,14 @@ class DiffusionOperator(TimeDependentOperator, LaplacianOperatorBase):
 
         assert self.dimensions == discr.dimensions
 
-        from grudge.mesh import check_bc_coverage
+        #from grudge.mesh import check_bc_coverage
         check_bc_coverage(discr.mesh, [self.dirichlet_tag, self.neumann_tag])
 
         return BoundDiffusionOperator(self, discr)
 
+    def m_sym_operator(self):
+        return self.sym_operator(apply_minv=True)
+
     def estimate_timestep(self, discr, 
             stepper=None, stepper_class=None, stepper_args=None,
             t=None, fields=None):
@@ -85,7 +87,8 @@ class DiffusionOperator(TimeDependentOperator, LaplacianOperatorBase):
 
 
 
-class BoundDiffusionOperator(grudge.iterative.OperatorBase):
+#class BoundDiffusionOperator(grudge.iterative.OperatorBase):
+class BoundDiffusionOperator():
     """Returned by :meth:`DiffusionOperator.bind`."""
 
     def __init__(self, diffusion_op, discr):
@@ -105,6 +108,8 @@ class BoundDiffusionOperator(grudge.iterative.OperatorBase):
         self.poincare_mean_value_hack = (
                 len(self.discr.get_boundary(BTAG_ALL).nodes)
                 == len(self.discr.get_boundary(diffusion_op.neumann_tag).nodes))
+            
+
 
     def __call__(self, t, u):
         dop = self.diffusion_op
diff --git a/grudge/models/euler2.py b/grudge/models/euler2.py
new file mode 100644
index 0000000000000000000000000000000000000000..5fb2663d86a0a96d5520a7338c764332815c9b38
--- /dev/null
+++ b/grudge/models/euler2.py
@@ -0,0 +1,211 @@
+# -*- coding: utf8 -*-
+"""Operators modeling advective phenomena."""
+
+from __future__ import division, absolute_import
+
+__copyright__ = "Copyright (C) 2009-2017 Andreas Kloeckner, Bogdan Enache"
+
+__license__ = """
+Permission is hereby granted, free of charge, to any person obtaining a copy
+of this software and associated documentation files (the "Software"), to deal
+in the Software without restriction, including without limitation the rights
+to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+copies of the Software, and to permit persons to whom the Software is
+furnished to do so, subject to the following conditions:
+
+The above copyright notice and this permission notice shall be included in
+all copies or substantial portions of the Software.
+
+THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
+THE SOFTWARE.
+"""
+
+
+import numpy as np
+import numpy
+import numpy.linalg as la
+
+from grudge.models import HyperbolicOperator
+from grudge import sym
+from pytools.obj_array import make_obj_array, join_fields
+
+
+
+class EquationOfState(object):
+    def q_to_p(self, op, q):
+        raise NotImplementedError
+
+    def p_to_e(self, p, rho, u):
+        raise NotImplementedError
+
+class GammaLawEOS(EquationOfState):
+    # FIXME Shouldn't gamma only occur in the equation of state?
+    # I.e. shouldn't all uses of gamma go through the EOS?
+
+    def __init__(self, gamma):
+        self.gamma = gamma
+
+    def q_to_p(self, op, q):
+        return (self.gamma-1)*(op.e(q)-0.5*numpy.dot(op.rho_u(q), op.u(q)))
+
+    def p_to_e(self, p, rho, u):
+        return p / (self.gamma - 1) + rho / 2 * numpy.dot(u, u)
+
+class PolytropeEOS(GammaLawEOS):
+    # inverse is same as superclass
+
+    def q_to_p(self, op, q):
+        return  op.rho(q)**self.gamma
+
+
+
+
+
+class EulerOperator(HyperbolicOperator):
+    def __init__(self, dimensions,
+            gamma=None, mu=0, bc_inflow=None,
+            bc_outflow=None,
+            bc_noslip=None,
+            bc_supersonic_inflow=None,
+            prandtl=None, spec_gas_const=1.0,
+            equation_of_state=None,
+            inflow_tag="inflow",
+            outflow_tag="outflow",
+            noslip_tag="noslip",
+            wall_tag="wall",
+            supersonic_inflow_tag="supersonic_inflow",
+            supersonic_outflow_tag="supersonic_outflow",
+            source=None,
+            second_order_scheme=None,
+            artificial_viscosity_mode=None,
+            quad_tag = "product"
+            ):
+        self.dimensions = dimensions
+        self.equation_of_state = equation_of_state
+        self.bc_inflow = sym.interp("vol",sym.BTAG_ALL)(bc_inflow)
+
+
+    def rho(self, q):
+        return q[0]
+
+    def e(self, q):
+        return q[1]
+
+    def rho_u(self, q):
+        return q[2:2+self.dimensions]
+
+    def u(self, q):
+        rho_u = self.rho_u(q)
+        return make_obj_array([
+                rho_u[i]/self.rho(q)
+                for i in range(self.dimensions)])
+
+    def p(self,q):
+        return self.equation_of_state.q_to_p(self, q)
+
+    def flux(self, q):
+        from pytools import delta
+        return [ # one entry for each flux direction
+                join_fields(
+                    # flux rho
+                    self.rho(q),
+
+                    # flux E
+                    0,
+
+                    # flux rho_u
+                    0, 0
+                    )
+                for i in range(self.dimensions)]
+
+        return [ # one entry for each flux direction
+                join_fields(
+                    # flux rho
+                    self.rho_u(q)[i],
+
+                    # flux E
+                    (self.e(q)+self.p(q))*self.u(q)[i],
+
+                    # flux rho_u
+                    make_obj_array([
+                        self.rho_u(q)[i]*self.u(q)[j]
+                        + delta(i,j) * self.p(q)
+                        for j in range(self.dimensions)
+                        ])
+                    )
+                for i in range(self.dimensions)]
+
+
+    def lax_flux(self, speed, state, flux):
+
+        normal = sym.normal(state.dd, self.dimensions)
+
+        wave_speed_ph = speed
+        state_ph = state
+        fluxes_ph = flux
+
+        from pymbolic.primitives import Max
+        #penalty = Max((wave_speed_ph.int,wave_speed_ph.ext))*(normal.ext*state_ph.ext+normal.int*state_ph.int)
+        penalty = 0
+
+        return state.avg*np.dot([1,1],normal)
+        #dotty = np.dot(normal, flux)
+        return (sum(n_i*(sym.int_tpair(f_i).avg) for n_i, f_i in zip(normal, fluxes_ph)))
+
+        #num_flux = 0.5*(sum(n_i*(sym.int_tpair(f_i).int+sym.int_tpair(f_i).ext) for n_i, f_i in zip(normal, fluxes_ph))
+        num_flux = (sum(n_i*(sym.int_tpair(f_i).avg) for n_i, f_i in zip(normal, fluxes_ph))
+            + penalty)
+        return num_flux
+
+    def characteristic_velocity(self,state):
+        from grudge.symbolic.operators import ElementwiseMaxOperator
+        from grudge.symbolic.primitives import CFunction
+        sqrt = CFunction("sqrt")
+
+        sound_speed = sqrt(self.equation_of_state.gamma*self.p(state)/self.rho(state))
+        u = self.u(state)
+        speed = sqrt(numpy.dot(u, u)) + sound_speed
+        return ElementwiseMaxOperator("vol","vol")(speed)
+
+
+        
+    def sym_operator(self):
+        q = sym.make_sym_array("q", self.dimensions+2)
+        def mflux(pair):
+            normal = sym.normal(pair.dd, self.dimensions)
+            return pair.avg*np.dot(normal, [1,1])
+
+        def flux(pair):
+            return sym.interp(pair.dd, "all_faces")(
+                    mflux(pair))
+        stiff_t = sym.stiffness_t(self.dimensions)
+        return  sym.InverseMassOperator()(
+                np.dot(stiff_t, [q,q])
+                - sym.FaceMassOperator()(
+                    flux(sym.int_tpair(q))
+                    + flux(sym.bv_tpair(sym.BTAG_ALL, q,self.bc_inflow ))
+                ))
+
+
+
+        def to_all_faces(pair, dd):
+            return sym.interp(dd, "all_faces")(pair)
+
+        speed = 0 #self.characteristic_velocity(q)
+        volq_flux = self.flux(q)
+        stiff_t = sym.stiffness_t(self.dimensions)
+        ans = sym.InverseMassOperator()(
+                np.dot(stiff_t, [q,q])
+                - sym.FaceMassOperator()(
+                    to_all_faces(self.lax_flux(sym.int_tpair(speed), sym.int_tpair(q), self.flux(q) ), sym.int_tpair(q).dd)
+                    #+ flux(sym.bv_tpair(sym.BTAG_ALL, q, bc_in))
+
+                    )
+                )
+        return join_fields(ans,speed)
+# vim: foldmethod=marker
diff --git a/grudge/models/gas_dynamics/__init__.py b/grudge/models/gas_dynamics/__init__.py
index e5a8ddc9418f1bac4776f252d2951fab3d3c1a2b..4ff7e83a13b3273a8c8b9c15ab2f49a482f23331 100644
--- a/grudge/models/gas_dynamics/__init__.py
+++ b/grudge/models/gas_dynamics/__init__.py
@@ -32,36 +32,19 @@ THE SOFTWARE.
 
 import numpy
 import grudge.tools
-import grudge.mesh
-import grudge.data
+import meshmode.mesh
 from grudge.models import TimeDependentOperator
 from pytools import Record
 from grudge.tools import is_zero
-from grudge.second_order import (
-        StabilizedCentralSecondDerivative,
-        CentralSecondDerivative,
-        IPDGSecondDerivative)
-from grudge.symbolic.primitives import make_common_subexpression as cse
+#from pymbolic.primitives import make_common_subexpression as cse
 from pytools import memoize_method
-from grudge.symbolic.tools import make_sym_vector
 from pytools.obj_array import make_obj_array, join_fields
+from grudge import sym
 
 AXES = ["x", "y", "z", "w"]
 
-from grudge.symbolic.operators import (
-        QuadratureGridUpsampler,
-        QuadratureInteriorFacesGridUpsampler)
-
-to_vol_quad = QuadratureGridUpsampler("gasdyn_vol")
-
-# It is recommended (though not required) that these two
-# remain the same so that they can be computed together
-# by the CUDA backend
-
-to_int_face_quad = QuadratureInteriorFacesGridUpsampler("gasdyn_face")
-to_bdry_quad = QuadratureGridUpsampler("gasdyn_face")
-
-
+def cse(x, stsr=None):
+    return x
 
 
 # {{{ equations of state
@@ -94,9 +77,6 @@ class PolytropeEOS(GammaLawEOS):
 # }}}
 
 
-
-
-
 class GasDynamicsOperator(TimeDependentOperator):
     """An nD Navier-Stokes and Euler operator.
 
@@ -141,8 +121,9 @@ class GasDynamicsOperator(TimeDependentOperator):
             supersonic_inflow_tag="supersonic_inflow",
             supersonic_outflow_tag="supersonic_outflow",
             source=None,
-            second_order_scheme=CentralSecondDerivative(),
+            second_order_scheme=None,
             artificial_viscosity_mode=None,
+            quad_tag = "product"
             ):
         """
         :param source: should implement
@@ -151,9 +132,6 @@ class GasDynamicsOperator(TimeDependentOperator):
 
         :param artificial_viscosity_mode:
         """
-        from grudge.data import (
-                TimeConstantGivenFunction,
-                ConstantGivenFunction)
 
         if gamma is not None:
             if equation_of_state is not None:
@@ -165,9 +143,10 @@ class GasDynamicsOperator(TimeDependentOperator):
 
             equation_of_state = GammaLawEOS(gamma)
 
-        dull_bc = TimeConstantGivenFunction(
-                ConstantGivenFunction(make_obj_array(
-                    [1, 1] + [0]*dimensions)))
+        mostly_one = sym.nodes(dimensions)/sym.nodes(dimensions)
+        mostly_zero = sym.nodes(dimensions)*0
+        #dull_bc = make_obj_array( [1, 1] + [0]*dimensions)
+        dull_bc = make_obj_array( [mostly_one,mostly_one] + [mostly_one]*dimensions)
         if bc_inflow is None:
             bc_inflow = dull_bc
         if bc_outflow is None:
@@ -183,6 +162,8 @@ class GasDynamicsOperator(TimeDependentOperator):
         self.spec_gas_const = spec_gas_const
         self.mu = mu
 
+        self.BAD_DEBUG = None
+
         self.bc_inflow = bc_inflow
         self.bc_outflow = bc_outflow
         self.bc_noslip = bc_noslip
@@ -205,6 +186,8 @@ class GasDynamicsOperator(TimeDependentOperator):
             raise ValueError("artificial_viscosity_mode has an invalid value")
 
         self.artificial_viscosity_mode = artificial_viscosity_mode
+        self.quad_tag = None #quad_tag
+
 
 
 
@@ -212,14 +195,20 @@ class GasDynamicsOperator(TimeDependentOperator):
 
     # {{{ conversions ---------------------------------------------------------
     def state(self):
-        return make_sym_vector("q", self.dimensions+2)
+        return sym.make_sym_array("q", self.dimensions+2)
 
     @memoize_method
     def volq_state(self):
+        return self.state()
+        quad_dd = sym.DOFDesc("vol", self.quad_tag)
+        to_vol_quad = sym.interp("vol", quad_dd)
         return cse(to_vol_quad(self.state()), "vol_quad_state")
 
     @memoize_method
     def faceq_state(self):
+        int_faces_dd = sym.DOFDesc(sym.FACE_RESTR_INTERIOR, self.quad_tag)
+        to_int_face_quad = sym.interp("vol", int_faces_dd)
+        #return cse(sym.int_tpair(self.state(), self.quad_tag), "face_quad_state")
         return cse(to_int_face_quad(self.state()), "face_quad_state")
 
     @memoize_method
@@ -256,12 +245,6 @@ class GasDynamicsOperator(TimeDependentOperator):
     def cse_p(self, q):
         return cse(self.p(q), "p")
 
-    def temperature(self, q):
-        c_v = 1 / (self.equation_of_state.gamma - 1) * self.spec_gas_const
-        return (self.e(q)/self.rho(q) - 0.5 * numpy.dot(self.u(q), self.u(q))) / c_v
-
-    def cse_temperature(self, q):
-        return cse(self.temperature(q), "temperature")
 
     def get_mu(self, q, to_quad_op):
         """
@@ -298,7 +281,7 @@ class GasDynamicsOperator(TimeDependentOperator):
 
     def primitive_to_conservative(self, prims, use_cses=True):
         if not use_cses:
-            from grudge.symbolic.primitives import make_common_subexpression as cse
+            def cse(x, name): return sym.cse(x,name)
         else:
             def cse(x, name): return x
 
@@ -314,7 +297,7 @@ class GasDynamicsOperator(TimeDependentOperator):
 
     def conservative_to_primitive(self, q, use_cses=True):
         if use_cses:
-            from grudge.symbolic.primitives import make_common_subexpression as cse
+            def cse(x, name): return sym.cse(x,name)
         else:
             def cse(x, name): return x
 
@@ -332,12 +315,14 @@ class GasDynamicsOperator(TimeDependentOperator):
         sound_speed = cse(sqrt(
             self.equation_of_state.gamma*self.cse_p(state)/self.cse_rho(state)),
             "sound_speed")
+        #sound_speed = sqrt(self.cse_p(state))
+        #sound_speed = self.cse_p(state)
         u = self.cse_u(state)
         speed = cse(sqrt(numpy.dot(u, u)), "norm_u") + sound_speed
-        return ElementwiseMaxOperator()(speed)
+        return ElementwiseMaxOperator("vol","vol")(speed)
 
     def bind_characteristic_velocity(self, discr):
-        state = make_sym_vector("q", self.dimensions+2)
+        state = sym.make_sym_array("q", self.dimensions+2)
 
         compiled = discr.compile(
                 self.characteristic_velocity_optemplate(state))
@@ -347,116 +332,7 @@ class GasDynamicsOperator(TimeDependentOperator):
 
         return do
 
-    # }}}
-
-    # {{{ helpers for second-order part ---------------------------------------
-
-    # {{{ compute gradient of state ---------------------------------------
-    def grad_of(self, var, faceq_var):
-        from grudge.second_order import SecondDerivativeTarget
-        grad_tgt = SecondDerivativeTarget(
-                self.dimensions, strong_form=False,
-                operand=var,
-                int_flux_operand=faceq_var,
-                bdry_flux_int_operand=faceq_var)
-
-        self.second_order_scheme.grad(grad_tgt,
-                bc_getter=self.get_boundary_condition_for,
-                dirichlet_tags=self.get_boundary_tags(),
-                neumann_tags=[])
-
-        return grad_tgt.minv_all
-
-    def grad_of_state(self):
-        dimensions = self.dimensions
-
-        state = self.state()
-
-        dq = numpy.zeros((len(state), dimensions), dtype=object)
-
-        for i in range(len(state)):
-            dq[i,:] = self.grad_of(
-                    state[i], self.faceq_state()[i])
-
-        return dq
-
-    def grad_of_state_func(self, func, of_what_descr):
-        return cse(self.grad_of(
-            func(self.volq_state()),
-            func(self.faceq_state())),
-            "grad_"+of_what_descr)
-
-    # }}}
-
-    # {{{ viscous stress tensor
-
-    def tau(self, to_quad_op, state, mu=None):
-        faceq_state = self.faceq_state()
-
-        dimensions = self.dimensions
-
-        # {{{ compute gradient of u ---------------------------------------
-        # Use the product rule to compute the gradient of
-        # u from the gradient of (rho u). This ensures we don't
-        # compute the derivatives twice.
-
-        from pytools.obj_array import with_object_array_or_scalar
-        dq = with_object_array_or_scalar(
-                to_quad_op, self.grad_of_state())
-
-        q = cse(to_quad_op(state))
-
-        du = numpy.zeros((dimensions, dimensions), dtype=object)
-        for i in range(dimensions):
-            for j in range(dimensions):
-                du[i,j] = cse(
-                        (dq[i+2,j] - self.cse_u(q)[i] * dq[0,j]) / self.rho(q),
-                        "du%d_d%s" % (i, AXES[j]))
-
-        # }}}
-
-        # {{{ put together viscous stress tau -----------------------------
-        from pytools import delta
-
-        if mu is None:
-            mu = self.get_mu(q, to_quad_op)
-
-        tau = numpy.zeros((dimensions, dimensions), dtype=object)
-        for i in range(dimensions):
-            for j in range(dimensions):
-                tau[i,j] = cse(mu * cse(du[i,j] + du[j,i] -
-                           2/self.dimensions * delta(i,j) * numpy.trace(du)),
-                           "tau_%d%d" % (i, j))
-
-        return tau
-
-        # }}}
-
-    # }}}
-
-    # }}}
-
-    # {{{ heat conduction
-
-    def heat_conduction_coefficient(self, to_quad_op):
-        mu = self.get_mu(self.state(), to_quad_op)
-        if self.prandtl is None or numpy.isinf(self.prandtl):
-            return 0
-
-        eos = self.equation_of_state
-        return (mu / self.prandtl) * (eos.gamma / (eos.gamma-1))
-
-    def heat_conduction_grad(self, to_quad_op):
-        grad_p_over_rho = self.grad_of_state_func(
-                lambda state: self.p(state)/self.rho(state),
-                "p_over_rho")
-
-        return (self.heat_conduction_coefficient(to_quad_op)
-                * to_quad_op(grad_p_over_rho))
-
-    # }}}
 
-    # {{{ flux
 
     def flux(self, q):
         from pytools import delta
@@ -488,9 +364,15 @@ class GasDynamicsOperator(TimeDependentOperator):
           BC is linearized.
         """
         if state0 is None:
-            state0 = make_sym_vector(bc_name, self.dimensions+2)
+            raise ValueError
+            #state0 = sym.make_sym_array(bc_name, self.dimensions+2)
 
-        state0 = cse(to_bdry_quad(state0))
+        #boundary_dd = sym.DOFDesc(sym.BTAG_ALL, self.quad_tag)
+        boundary_dd = sym.DOFDesc(sym.DTAG_BOUNDARY(tag), self.quad_tag)
+        to_bdry_quad = sym.interp(sym.BTAG_ALL, boundary_dd)
+        #state0 = cse(to_bdry_quad(state0))
+        to_bdry_quad_f_a = sym.interp("vol", boundary_dd)
+        state0 = cse(to_bdry_quad_f_a(state0))
 
         rho0 = self.rho(state0)
         p0 = self.cse_p(state0)
@@ -498,8 +380,8 @@ class GasDynamicsOperator(TimeDependentOperator):
 
         c0 = (self.equation_of_state.gamma * p0 / rho0)**0.5
 
-        from grudge.symbolic import RestrictToBoundary
-        bdrize_op = RestrictToBoundary(tag)
+        #bdrize_op = sym.interp("vol",sym.BTAG_ALL)
+        bdrize_op = sym.interp("vol",boundary_dd)
 
         class SingleBCInfo(Record):
             pass
@@ -516,9 +398,10 @@ class GasDynamicsOperator(TimeDependentOperator):
                 - p0, "dpm"))
 
     def outflow_state(self, state):
-        from grudge.symbolic import make_normal
-        normal = make_normal(self.outflow_tag, self.dimensions)
-        bc = self.make_bc_info("bc_q_out", self.outflow_tag, state)
+
+        normal = sym.normal(sym.DTAG_BOUNDARY(self.outflow_tag), self.dimensions)
+        #bc = self.make_bc_info("bc_q_out", self.outflow_tag, state, self.bc_outflow)
+        bc = self.make_bc_info("bc_q_out", self.outflow_tag, state, state)
 
         # see grudge/doc/maxima/euler.mac
         return join_fields(
@@ -552,29 +435,24 @@ class GasDynamicsOperator(TimeDependentOperator):
             + normal*numpy.dot(normal, bc.dumvec)/2 + bc.dpm*normal/(2*bc.c0*bc.rho0), "bc_u_"+name))
 
     def inflow_state(self, state):
-        from grudge.symbolic import make_normal
-        normal = make_normal(self.inflow_tag, self.dimensions)
-        bc = self.make_bc_info("bc_q_in", self.inflow_tag, state)
+        normal = sym.normal(sym.DTAG_BOUNDARY(self.inflow_tag), self.dimensions)
+        #bc = self.make_bc_info("bc_q_in", self.inflow_tag, state, self.bc_inflow)
+        bc = self.make_bc_info("bc_q_in", self.inflow_tag, state, state)
         return self.inflow_state_inner(normal, bc, "inflow")
 
     def noslip_state(self, state):
-        from grudge.symbolic import make_normal
-        state0 = join_fields(
-            make_sym_vector("bc_q_noslip", 2),
-            [0]*self.dimensions)
-        normal = make_normal(self.noslip_tag, self.dimensions)
-        bc = self.make_bc_info("bc_q_noslip", self.noslip_tag, state, state0)
+        normal = sym.normal(sym.DTAG_BOUNDARY(self.noslip_tag), self.dimensions)
+        #bc = self.make_bc_info("bc_q_noslip", self.noslip_tag, state, self.bc_noslip)
+        bc = self.make_bc_info("bc_q_noslip", self.noslip_tag, state, state)
         return self.inflow_state_inner(normal, bc, "noslip")
 
     def wall_state(self, state):
-        from grudge.symbolic import RestrictToBoundary
-        bc = RestrictToBoundary(self.wall_tag)(state)
+        bc = sym.interp("vol",sym.BTAG_ALL)(state)
         wall_rho = self.rho(bc)
         wall_e = self.e(bc) # <3 eve
         wall_rho_u = self.rho_u(bc)
 
-        from grudge.symbolic import make_normal
-        normal = make_normal(self.wall_tag, self.dimensions)
+        normal = sym.normal(sym.DTAG_BOUNDARY(self.wall_tag), self.dimensions)
 
         return join_fields(
                 wall_rho,
@@ -596,18 +474,19 @@ class GasDynamicsOperator(TimeDependentOperator):
     def get_conservative_boundary_conditions(self):
         state = self.state()
 
-        from grudge.symbolic import RestrictToBoundary
         return {
-                self.supersonic_inflow_tag:
-                make_sym_vector("bc_q_supersonic_in", self.dimensions+2),
+                self.supersonic_inflow_tag: self.bc_supersonic_inflow,
+                #sym.make_sym_array("bc_q_supersonic_in", self.dimensions+2),
                 self.supersonic_outflow_tag:
-                RestrictToBoundary(self.supersonic_outflow_tag)(
-                            (state)),
+                #sym.interp("vol",self.supersonic_outflow_tag)(state), #TODO: It should be this according to my em translation <- but (it spits out an error)
+                sym.interp("vol",sym.DTAG_BOUNDARY(self.supersonic_outflow_tag))(state), #TODO: Attempt 2
+                #sym.interp("vol",sym.BTAG_ALL)(state),
                 self.wall_tag: self.wall_state(state),
                 }
 
     @memoize_method
     def get_boundary_tags(self):
+        return set([self.inflow_tag,])
         return (set(self.get_primitive_boundary_conditions().keys())
                 | set(self.get_conservative_boundary_conditions().keys()))
 
@@ -615,11 +494,10 @@ class GasDynamicsOperator(TimeDependentOperator):
     def _normalize_expr(self, expr):
         """Normalize expressions for use as hash keys."""
         from grudge.symbolic.mappers import (
-                QuadratureUpsamplerRemover,
+                QuadratureRemover,
                 CSERemover)
 
-        return CSERemover()(
-                QuadratureUpsamplerRemover({}, do_warn=False)(expr))
+        return CSERemover()(QuadratureRemover()(expr))
 
     @memoize_method
     def _get_norm_primitive_exprs(self):
@@ -632,9 +510,11 @@ class GasDynamicsOperator(TimeDependentOperator):
     def get_boundary_condition_for(self, tag, expr):
         prim_bcs = self.get_primitive_boundary_conditions()
         cons_bcs = self.get_conservative_boundary_conditions()
+        boundary_dd = sym.DOFDesc(sym.BTAG_ALL, self.quad_tag)
+        to_bdry_quad = sym.interp(sym.DD_VOLUME, boundary_dd)
 
         if tag in prim_bcs:
-            # BC is given in primitive variables, avoid converting
+            # BC is giVEn in primitive variables, avoid converting
             # to conservative and back.
             try:
                 norm_expr = self._normalize_expr(expr)
@@ -648,74 +528,71 @@ class GasDynamicsOperator(TimeDependentOperator):
             # BC is given in conservative variables, no potential
             # for optimization.
 
-            cbstate = to_bdry_quad(cons_bcs[tag])
+            #cbstate = to_bdry_quad(cons_bcs[tag])
+            cbstate = cons_bcs[tag]
 
         # 'cbstate' is the boundary state in conservative variables.
 
-        from grudge.symbolic.mappers import QuadratureUpsamplerRemover
-        expr = QuadratureUpsamplerRemover({}, do_warn=False)(expr)
+        from grudge.symbolic.mappers import QuadratureRemover
+        expr = QuadratureRemover()(expr)
 
         def subst_func(expr):
             from pymbolic.primitives import Subscript, Variable
+            #boundary_dd = sym.DOFDesc(sym.BTAG_ALL, self.quad_tag)
+            #to_bdry_quad = sym.interp(sym.DD_VOLUME, boundary_dd)
 
             if isinstance(expr, Subscript):
                 assert (isinstance(expr.aggregate, Variable) 
                         and expr.aggregate.name == "q")
 
                 return cbstate[expr.index]
-            elif isinstance(expr, Variable) and expr.name =="sensor":
-                from grudge.symbolic import RestrictToBoundary
-                result = RestrictToBoundary(tag)(self.sensor())
-                return cse(to_bdry_quad(result), "bdry_sensor")
 
-        from grudge.symbolic import SubstitutionMapper
+        from grudge.symbolic.mappers import SubstitutionMapper
         return SubstitutionMapper(subst_func)(expr)
 
     # }}}
 
-    # {{{ second order part
-    def div(self, vol_operand, int_face_operand):
-        from grudge.second_order import SecondDerivativeTarget
-        div_tgt = SecondDerivativeTarget(
-                self.dimensions, strong_form=False,
-                operand=vol_operand,
-                int_flux_operand=int_face_operand)
 
-        self.second_order_scheme.div(div_tgt,
-                bc_getter=self.get_boundary_condition_for,
-                dirichlet_tags=list(self.get_boundary_tags()),
-                neumann_tags=[])
 
-        return div_tgt.minv_all
+    # {{{ operator template ---------------------------------------------------
 
-    def make_second_order_part(self):
-        state = self.state()
-        faceq_state = self.faceq_state()
-        volq_state = self.volq_state()
+    def make_lax_friedrichs_flux(self, wave_speed, state, fluxes, bdry_tags_states_and_fluxes,strong):
 
-        volq_tau_mat = self.tau(to_vol_quad, state)
-        faceq_tau_mat = self.tau(to_int_face_quad, state)
+        def lff(fvph):
+            n = len(state)
+            d = len(fluxes)
+            normal = sym.normal(fvph.dd, d)
 
-        return join_fields(
-                0, 
-                self.div(
-                    numpy.sum(volq_tau_mat*self.cse_u(volq_state), axis=1)
-                    + self.heat_conduction_grad(to_vol_quad)
-                    ,
-                    numpy.sum(faceq_tau_mat*self.cse_u(faceq_state), axis=1)
-                    + self.heat_conduction_grad(to_int_face_quad)
-                    ,
-                    ),
-                [
-                    self.div(volq_tau_mat[i], faceq_tau_mat[i])
-                    for i in range(self.dimensions)]
-                )
+            wave_speed_ph = fvph[0]
+            state_ph = fvph[1:1+n]
+            fluxes_ph = [fvph[1+i*n:1+(i+1)*n] for i in range(1, d+1)]
+            print(fvph.dd)
+
+            from pymbolic.primitives import Max
+            #penalty = Max((wave_speed_ph.int,wave_speed_ph.ext))*(state_ph.ext-state_ph.int)
+            penalty = 0
+
+            if not strong:
+                num_flux = 0.5*(sum(n_i*(f_i.int+f_i.ext) for n_i, f_i in zip(normal, fluxes_ph))
+                - penalty)
+            else:
+                num_flux = 0.5*(sum(n_i*(f_i.int-f_i.ext) for n_i, f_i in zip(normal, fluxes_ph))
+                    + penalty)
+            return num_flux
+
+        all_faces_dd = sym.DOFDesc(sym.FACE_RESTR_ALL, self.quad_tag)
+        def fluxy(pair):
+            return sym.interp(pair.dd, all_faces_dd)(lff(pair))
+
+        int_operand = join_fields(wave_speed, state, *fluxes)
+        
+
+        #return fluxy(sym.int_tpair(int_operand, self.quad_tag))
+        return (fluxy(sym.int_tpair(int_operand, self.quad_tag)) + sum(
+                fluxy(sym.bv_tpair(sym.DOFDesc(sym.DTAG_BOUNDARY(tag), self.quad_tag),int_operand,join_fields(1, bdry_state, *bdry_fluxes) ))
+            for tag, bdry_state, bdry_fluxes in bdry_tags_states_and_fluxes))
 
-    # }}}
 
-    # {{{ operator template ---------------------------------------------------
-    def make_extra_terms(self):
-        return 0
 
     def sym_operator(self, sensor_scaling=None, viscosity_only=False):
         u = self.cse_u
@@ -723,22 +600,9 @@ class GasDynamicsOperator(TimeDependentOperator):
         rho_u = self.rho_u
         p = self.p
         e = self.e
+        int_faces_dd = sym.DOFDesc(sym.FACE_RESTR_INTERIOR, self.quad_tag)
+        to_int_face_quad = sym.interp("vol", int_faces_dd)
 
-        # {{{ artificial diffusion
-        def make_artificial_diffusion():
-            if self.artificial_viscosity_mode not in ["diffusion"]:
-                return 0
-
-            dq = self.grad_of_state()
-
-            return make_obj_array([
-                self.div(
-                    to_vol_quad(self.sensor())*to_vol_quad(dq[i]),
-                    to_int_face_quad(self.sensor())*to_int_face_quad(dq[i])) 
-                for i in range(dq.shape[0])])
-        # }}}
-
-        # {{{ state setup
 
         volq_flux = self.flux(self.volq_state())
         faceq_flux = self.flux(self.faceq_state())
@@ -751,13 +615,10 @@ class GasDynamicsOperator(TimeDependentOperator):
         has_viscosity = not is_zero(self.get_mu(self.state(), to_quad_op=None))
 
         # }}}
+        boundary_dd = sym.DOFDesc(sym.BTAG_ALL, self.quad_tag)
+        to_bdry_quad = sym.interp(sym.DD_VOLUME, boundary_dd)
 
         # {{{ operator assembly -----------------------------------------------
-        from grudge.flux.tools import make_lax_friedrichs_flux
-        from grudge.symbolic.operators import InverseMassOperator
-
-        from grudge.symbolic.tools import make_stiffness_t
-
         primitive_bcs_as_quad_conservative = dict(
                 (tag, self.primitive_to_conservative(to_bdry_quad(bc)))
                 for tag, bc in 
@@ -769,153 +630,38 @@ class GasDynamicsOperator(TimeDependentOperator):
                 self.get_boundary_condition_for(tag, s_i) for s_i in state])
             return tag, bc, self.flux(bc)
 
-        first_order_part = InverseMassOperator()(
-                numpy.dot(make_stiffness_t(self.dimensions), volq_flux)
-                - make_lax_friedrichs_flux(
-                    wave_speed=cse(to_int_face_quad(speed), "emax_c"),
 
-                    state=self.faceq_state(), fluxes=faceq_flux,
-                    bdry_tags_states_and_fluxes=[
-                        get_bc_tuple(tag) for tag in self.get_boundary_tags()],
-                    strong=False))
+        quad_dd = sym.DOFDesc("vol", self.quad_tag)
+        all_faces_dd = sym.DOFDesc(sym.FACE_RESTR_ALL, self.quad_tag)
+        stiff_t = sym.stiffness_t(self.dimensions, quad_dd, "vol")
+
+        first_order_part = sym.InverseMassOperator()(
+                numpy.dot(stiff_t, volq_flux)
+                - sym.FaceMassOperator(all_faces_dd, "vol")(self.make_lax_friedrichs_flux(
+                    speed,
+                    self.state(),
+                    self.flux(self.state()), #faceq_flux,
+                    [get_bc_tuple(tag) for tag in self.get_boundary_tags()],
+                    False))
+                )
 
         if viscosity_only:
             first_order_part = 0*first_order_part
 
         result = join_fields(
-                first_order_part 
-                + self.make_second_order_part()
-                + make_artificial_diffusion()
-                + self.make_extra_terms(),
+                first_order_part ,
+                #+ self.make_second_order_part()
+                #+ make_artificial_diffusion()
+                #+ self.make_extra_terms(),
                  speed)
 
         if self.source is not None:
             result = result + join_fields(
-                    make_sym_vector("source_vect", len(self.state())),
+                    sym.make_sym_array("source_vect", len(self.state())),
                     # extra field for speed
                     0)
 
         return result
 
-        # }}}
-
-    # }}}
-
-    # {{{ operator binding ----------------------------------------------------
-    def bind(self, discr, sensor=None, sensor_scaling=None, viscosity_only=False):
-        if (sensor is None and 
-                self.artificial_viscosity_mode is not None):
-            raise ValueError("must specify a sensor if using "
-                    "artificial viscosity")
-
-        bound_op = discr.compile(self.sym_operator(
-            sensor_scaling=sensor_scaling,
-            viscosity_only=False))
-
-        from grudge.mesh import check_bc_coverage
-        check_bc_coverage(discr.mesh, [
-            self.inflow_tag,
-            self.outflow_tag,
-            self.noslip_tag,
-            self.wall_tag,
-            self.supersonic_inflow_tag,
-            self.supersonic_outflow_tag,
-            ])
-
-        if self.mu == 0 and not discr.get_boundary(self.noslip_tag).is_empty():
-            raise RuntimeError("no-slip BCs only make sense for "
-                    "viscous problems")
-
-        def rhs(t, q):
-            extra_kwargs = {}
-            if self.source is not None:
-                extra_kwargs["source_vect"] = self.source.volume_interpolant(
-                        t, q, discr)
-
-            if sensor is not None:
-                extra_kwargs["sensor"] = sensor(q)
-
-            opt_result = bound_op(q=q,
-                    bc_q_in=self.bc_inflow.boundary_interpolant(
-                        t, discr, self.inflow_tag),
-                    bc_q_out=self.bc_outflow.boundary_interpolant(
-                        t, discr, self.outflow_tag),
-                    bc_q_noslip=self.bc_noslip.boundary_interpolant(
-                        t, discr, self.noslip_tag),
-                    bc_q_supersonic_in=self.bc_supersonic_inflow
-                    .boundary_interpolant(t, discr, 
-                        self.supersonic_inflow_tag),
-                    **extra_kwargs
-                    )
-
-            max_speed = opt_result[-1]
-            ode_rhs = opt_result[:-1]
-            return ode_rhs, discr.nodewise_max(max_speed)
-
-        return rhs
-
-    # }}}
-
-    # {{{ timestep estimation -------------------------------------------------
-
-    def estimate_timestep(self, discr, 
-            stepper=None, stepper_class=None, stepper_args=None,
-            t=None, max_eigenvalue=None):
-        u"""Estimate the largest stable timestep, given a time stepper
-        `stepper_class`. If none is given, RK4 is assumed.
-        """
-
-        dg_factor = (discr.dt_non_geometric_factor()
-                * discr.dt_geometric_factor())
-
-        # see JSH/TW, eq. (7.32)
-        rk4_dt = dg_factor / (max_eigenvalue + self.mu / dg_factor)
-
-        from grudge.timestep.stability import \
-                approximate_rk4_relative_imag_stability_region
-        return rk4_dt * approximate_rk4_relative_imag_stability_region(
-                stepper, stepper_class, stepper_args)
-
-    # }}}
-
-
-
-
-# {{{ limiter (unfinished, deprecated)
-class SlopeLimiter1NEuler:
-    def __init__(self, discr, gamma, dimensions, op):
-        """Construct a limiter from Jan's book page 225
-        """
-        self.discr = discr
-        self.gamma=gamma
-        self.dimensions=dimensions
-        self.op=op
-
-        from grudge.symbolic.operators import AveragingOperator
-        self.get_average = AveragingOperator().bind(discr)
-
-    def __call__(self, fields):
-        from grudge.tools import join_fields
-
-        #get conserved fields
-        rho=self.op.rho(fields)
-        e=self.op.e(fields)
-        rho_velocity=self.op.rho_u(fields)
-
-        #get primitive fields
-        #to do
-
-        #reset field values to cell average
-        rhoLim=self.get_average(rho)
-        eLim=self.get_average(e)
-        temp=join_fields([self.get_average(rho_vel)
-                for rho_vel in rho_velocity])
-
-        #should do for primitive fields too
-
-        return join_fields(rhoLim, eLim, temp)
-
-# }}}
-
-
+    
 # vim: foldmethod=marker
diff --git a/grudge/models/poisson.py b/grudge/models/poisson.py
index a9f7145744c2c1a5e6e927680f341c3dbbf3e3e7..81e4b56abba9cbc34da6b37f3ce32d3b55d4560e 100644
--- a/grudge/models/poisson.py
+++ b/grudge/models/poisson.py
@@ -30,9 +30,9 @@ THE SOFTWARE.
 import numpy as np
 
 from grudge.models import Operator
-from grudge.second_order import LDGSecondDerivative
-import grudge.data
-import grudge.iterative
+from grudge.models.second_order import LDGSecondDerivative
+from grudge import sym
+#import grudge.iterative
 
 
 class LaplacianOperatorBase(object):
@@ -48,15 +48,14 @@ class LaplacianOperatorBase(object):
           :class:`grudge.models.diffusion.DiffusionOperator` needs this.
         """
 
-        from grudge.symbolic import Field, make_sym_vector
-        from grudge.second_order import SecondDerivativeTarget
+        from grudge.models.second_order import SecondDerivativeTarget
 
         if u is None:
-            u = Field("u")
+            u = sym.var("u")
         if dir_bc is None:
-            dir_bc = Field("dir_bc")
+            dir_bc = sym.var("dir_bc")
         if neu_bc is None:
-            neu_bc = Field("neu_bc")
+            neu_bc = sym.var("neu_bc")
 
         # strong_form here allows IPDG to reuse the value of grad u.
         grad_tgt = SecondDerivativeTarget(
@@ -66,10 +65,11 @@ class LaplacianOperatorBase(object):
         def grad_bc_getter(tag, expr):
             assert tag == self.dirichlet_tag
             return dir_bc
-        self.scheme.grad(grad_tgt,
+
+        self.scheme.grad(u, grad_tgt,
                 bc_getter=grad_bc_getter,
                 dirichlet_tags=[self.dirichlet_tag],
-                neumann_tags=[self.neumann_tag])
+                neumann_tags=[self.neumann_tag] )
 
         def apply_diff_tensor(v):
             if isinstance(self.diffusion_tensor, np.ndarray):
@@ -115,9 +115,9 @@ class PoissonOperator(Operator, LaplacianOperatorBase):
     """
 
     def __init__(self, dimensions, diffusion_tensor=None,
-            dirichlet_bc=grudge.data.ConstantGivenFunction(),
+            dirichlet_bc=0,
             dirichlet_tag="dirichlet",
-            neumann_bc=grudge.data.ConstantGivenFunction(),
+            neumann_bc=0,
             neumann_tag="neumann",
             scheme=LDGSecondDerivative()):
         self.dimensions = dimensions
@@ -145,11 +145,12 @@ class PoissonOperator(Operator, LaplacianOperatorBase):
         return BoundPoissonOperator(self, discr)
 
 
-class BoundPoissonOperator(grudge.iterative.OperatorBase):
+#class BoundPoissonOperator(grudge.iterative.OperatorBase):
+class BoundPoissonOperator():
     """Returned by :meth:`PoissonOperator.bind`."""
 
     def __init__(self, poisson_op, discr):
-        grudge.iterative.OperatorBase.__init__(self)
+        #grudge.iterative.OperatorBase.__init__(self)
         self.discr = discr
 
         pop = self.poisson_op = poisson_op
diff --git a/grudge/models/second_order.py b/grudge/models/second_order.py
index 24260c0f7012ab7bbd6f7adc0fc23a0248bbee32..b62cc9886d84a0c6ac1e6269b816699696012595 100644
--- a/grudge/models/second_order.py
+++ b/grudge/models/second_order.py
@@ -257,9 +257,12 @@ class SecondDerivativeTarget(object):
         else:
             return vec*operand
 
-    def normal_times_flux(self, flux):
-        from grudge.flux import make_normal
-        return self.vec_times(make_normal(self.dimensions), flux)
+    def normal_times_flux(self, flux, dd=None):
+        
+        if dd is None:
+            return self.vec_times(sym.normal(flux.dd, self.dimensions), flux)
+        else:
+            return self.vec_times(sym.normal(dd, self.dimensions), flux)
 
     def apply_diff(self, nabla, operand):
         from pytools.obj_array import make_obj_array, is_obj_array
@@ -275,11 +278,9 @@ class SecondDerivativeTarget(object):
 
     def _local_nabla(self):
         if self.strong_form:
-            from grudge.symbolic import make_stiffness
-            return make_stiffness(self.dimensions)
+            return sym.stiffness(self.dimensions)
         else:
-            from grudge.symbolic import make_stiffness_t
-            return make_stiffness_t(self.dimensions)
+            return sym.stiffness_t(self.dimensions)
 
     def add_derivative(self, operand=None):
         if operand is None:
@@ -288,10 +289,15 @@ class SecondDerivativeTarget(object):
         self.add_local_derivatives(
                 self.apply_diff(self._local_nabla(), operand))
 
-    def add_inner_fluxes(self, flux, expr):
-        from grudge.symbolic import get_flux_operator
+    def flux_face_interp(self, pair):
+        return sym.interp(pair.dd, "all_faces")(pair)
+
+    def add_inner_fluxes(self, flux):
+        #from grudge.symbolic import get_flux_operator
+        #self.inner_fluxes = self.inner_fluxes \
+                #+ get_flux_operator(self.strong_neg*flux)(expr)
         self.inner_fluxes = self.inner_fluxes \
-                + get_flux_operator(self.strong_neg*flux)(expr)
+                + self.flux_face_interp(self.strong_neg*flux)
 
     def add_boundary_flux(self, flux, volume_expr, bdry_expr, tag):
         from grudge.symbolic import BoundaryPair, get_flux_operator
@@ -309,10 +315,9 @@ class SecondDerivativeTarget(object):
 
     @property
     def minv_all(self):
-        from grudge.symbolic.primitives import make_common_subexpression as cse
         from grudge.symbolic.operators import InverseMassOperator
-        return (cse(InverseMassOperator()(self.local_derivatives), "grad_loc")
-                + cse(InverseMassOperator()(self.fluxes), "grad_flux"))
+        return (sym.cse(InverseMassOperator()(self.local_derivatives), "grad_loc")
+                + sym.cse(InverseMassOperator()(self.fluxes), "grad_flux"))
 
 # }}}
 
@@ -320,7 +325,7 @@ class SecondDerivativeTarget(object):
 # {{{ second derivative schemes
 
 class SecondDerivativeBase(object):
-    def grad(self, tgt, bc_getter, dirichlet_tags, neumann_tags):
+    def grad(self,u, tgt, bc_getter, dirichlet_tags, neumann_tags):
         """
         :param bc_getter: a function (tag, volume_expr) -> boundary expr.
           *volume_expr* will be None to query the Neumann condition.
@@ -328,19 +333,16 @@ class SecondDerivativeBase(object):
         n_times = tgt.normal_times_flux
 
         if tgt.strong_form:
-            def adjust_flux(f):
-                return n_times(u.int) - f
+            def adjust_flux(u, f,dd):
+                return n_times(u.int ,dd) - f
         else:
-            def adjust_flux(f):
+            def adjust_flux(u, f, dd):
                 return f
-
-        from grudge.flux import FluxScalarPlaceholder
-        u = FluxScalarPlaceholder()
+        int_u = sym.int_tpair(u)
 
         tgt.add_derivative()
         tgt.add_inner_fluxes(
-                adjust_flux(self.grad_interior_flux(tgt, u)),
-                tgt.int_flux_operand)
+                adjust_flux(int_u, self.grad_interior_flux(tgt, int_u), int_u.dd))
 
         for tag in dirichlet_tags:
             tgt.add_boundary_flux(
@@ -395,13 +397,12 @@ class LDGSecondDerivative(SecondDerivativeBase):
         return np.array([self.beta_value]*tgt.dimensions, dtype=np.float64)
 
     def grad_interior_flux(self, tgt, u):
-        from grudge.symbolic.primitives import make_common_subexpression as cse
         n_times = tgt.normal_times_flux
         v_times = tgt.vec_times
 
         return n_times(
-                cse(u.avg, "u_avg")
-                - v_times(self.beta(tgt), n_times(u.int-u.ext)))
+                sym.cse(u.avg, "u_avg")
+                - v_times(self.beta(tgt), n_times(u.int-u.ext, u.dd)), u.dd)
 
     def div(self, tgt, bc_getter, dirichlet_tags, neumann_tags):
         """
@@ -409,7 +410,6 @@ class LDGSecondDerivative(SecondDerivativeBase):
           *volume_expr* will be None to query the Neumann condition.
         """
 
-        from grudge.symbolic.primitives import make_common_subexpression as cse
         from grudge.flux import FluxVectorPlaceholder, PenaltyTerm
 
         n_times = tgt.normal_times_flux
@@ -431,13 +431,13 @@ class LDGSecondDerivative(SecondDerivativeBase):
 
         flux = n_times(flux_v.avg
                 + v_times(self.beta(tgt),
-                    cse(n_times(flux_v.int - flux_v.ext), "jump_v"))
+                    sym.cse(n_times(flux_v.int - flux_v.ext), "jump_v"))
                 - stab_term)
 
         from pytools.obj_array import make_obj_array
-        flux_arg_int = cse(make_obj_array(stab_term_generator.flux_args))
+        flux_arg_int = sym.cse(make_obj_array(stab_term_generator.flux_args))
 
-        tgt.add_derivative(cse(tgt.operand))
+        tgt.add_derivative(sym.cse(tgt.operand))
         tgt.add_inner_fluxes(adjust_flux(flux), flux_arg_int)
 
         self.add_div_bcs(tgt, bc_getter, dirichlet_tags, neumann_tags,
@@ -459,9 +459,8 @@ class IPDGSecondDerivative(SecondDerivativeBase):
         self.stab_coefficient = stab_coefficient
 
     def grad_interior_flux(self, tgt, u):
-        from grudge.symbolic.primitives import make_common_subexpression as cse
         n_times = tgt.normal_times_flux
-        return n_times(cse(u.avg, "u_avg"))
+        return n_times(sym.cse(u.avg, "u_avg"))
 
     def div(self, tgt, bc_getter, dirichlet_tags, neumann_tags):
         """
@@ -469,7 +468,6 @@ class IPDGSecondDerivative(SecondDerivativeBase):
           *volume_expr* will be None to query the Neumann condition.
         """
 
-        from grudge.symbolic.primitives import make_common_subexpression as cse
         from grudge.flux import FluxVectorPlaceholder, PenaltyTerm
 
         n_times = tgt.normal_times_flux
@@ -496,9 +494,9 @@ class IPDGSecondDerivative(SecondDerivativeBase):
         flux = n_times(pure_diff_v.avg - stab_term)
 
         from pytools.obj_array import make_obj_array
-        flux_arg_int = cse(make_obj_array(stab_term_generator.flux_args))
+        flux_arg_int = sym.cse(make_obj_array(stab_term_generator.flux_args))
 
-        tgt.add_derivative(cse(tgt.operand))
+        tgt.add_derivative(sym.cse(tgt.operand))
         tgt.add_inner_fluxes(adjust_flux(flux), flux_arg_int)
 
         self.add_div_bcs(tgt, bc_getter, dirichlet_tags, neumann_tags,
diff --git a/grudge/symbolic/compiler.py b/grudge/symbolic/compiler.py
index c555cea0a56e947b4a94e0cd1736e8fa2b97ce63..063db1b534e4943161e5e0c45f88cb6a828f2000 100644
--- a/grudge/symbolic/compiler.py
+++ b/grudge/symbolic/compiler.py
@@ -103,10 +103,11 @@ class LoopyKernelInstruction(Instruction):
     @memoize_method
     def get_dependencies(self):
         from pymbolic.primitives import Variable, Subscript
+        dep_mapper = _make_dep_mapper(include_subscripts=False)
         return set(
-                v
+                dep
                 for v in self.kernel_descriptor.input_mappings.values()
-                if isinstance(v, (Variable, Subscript)))
+                for dep in dep_mapper(v))
 
     def __str__(self):
         knl_str = "\n".join(
@@ -493,6 +494,12 @@ class Code(object):
                     if pre_assign_check is not None:
                         pre_assign_check(target, value)
 
+                    import pyopencl.array
+                    import numpy as np
+                    if isinstance(value,  pyopencl.array.Array):
+                        if np.isnan(value.get()).any():
+                            print(value.get())
+                            pu.db
                     context[target] = value
 
                 futures.extend(new_futures)
diff --git a/grudge/symbolic/dofdesc_inference.py b/grudge/symbolic/dofdesc_inference.py
index 832c6a0371764cceef497e67a4a9c27300a5825a..942b805b622b38663f8ba75403eaedc6c4190bb0 100644
--- a/grudge/symbolic/dofdesc_inference.py
+++ b/grudge/symbolic/dofdesc_inference.py
@@ -99,6 +99,8 @@ class DOFDescInferenceMapper(RecursiveMapper, CSECachingMapperMixin):
         self.name_to_dofdesc = name_to_dofdesc
 
     def infer_for_name(self, name):
+        #if(name=='expr_29'):
+            #pu.db
         try:
             return self.name_to_dofdesc[name]
         except KeyError:
@@ -157,6 +159,10 @@ class DOFDescInferenceMapper(RecursiveMapper, CSECachingMapperMixin):
     def map_nodal_sum(self, expr, enclosing_prec):
         return DOFDesc(DTAG_SCALAR)
 
+    def map_max(self, expr):
+        return self.map_multi_child(expr, expr.children)
+
+
     map_nodal_max = map_nodal_sum
     map_nodal_min = map_nodal_sum
 
diff --git a/grudge/symbolic/mappers/__init__.py b/grudge/symbolic/mappers/__init__.py
index ca11cbf0d92c96b6c9d00a6e035b7b840c2a85cd..59e4b0843cc003b421d7c09b379d3df88030fbbd 100644
--- a/grudge/symbolic/mappers/__init__.py
+++ b/grudge/symbolic/mappers/__init__.py
@@ -135,6 +135,7 @@ class OperatorReducerMixin(LocalOpReducerMixin, FluxOpReducerMixin):
     map_nodal_sum = _map_op_base
     map_nodal_min = _map_op_base
     map_nodal_max = _map_op_base
+    map_elementwise_max = _map_op_base
 
     map_stiffness = _map_op_base
     map_diff = _map_op_base
@@ -183,6 +184,7 @@ class IdentityMapperMixin(LocalOpReducerMixin, FluxOpReducerMixin):
     map_nodal_sum = map_elementwise_linear
     map_nodal_min = map_elementwise_linear
     map_nodal_max = map_elementwise_linear
+    map_elementwise_max = map_elementwise_linear
 
     map_stiffness = map_elementwise_linear
     map_diff = map_elementwise_linear
@@ -515,19 +517,23 @@ class GlobalToReferenceMapper(CSECachingMapperMixin, IdentityMapper):
             if with_jacobian:
                 rec_field = jac_tag * rec_field
             return sum(
+                ref_class(rst_axis, dd_in=dd_in)(
                     sym.inverse_metric_derivative(
                         rst_axis, expr.op.xyz_axis,
-                        ambient_dim=self.ambient_dim, dim=self.dim) *
-                    ref_class(rst_axis, dd_in=dd_in)(rec_field)
-                    for rst_axis in range(self.dim))
+                        ambient_dim=self.ambient_dim, dim=self.dim, dd=dd_in) * rec_field)
+                for rst_axis in range(self.dim))
 
         if isinstance(expr.op, op.MassOperator):
             return op.RefMassOperator(expr.op.dd_in, expr.op.dd_out)(
                     jac_in * self.rec(expr.field))
 
         elif isinstance(expr.op, op.InverseMassOperator):
-            return op.RefInverseMassOperator(expr.op.dd_in, expr.op.dd_out)(
-                1/jac_in * self.rec(expr.field))
+            # Chan, Hewett, and Warburton. "Weight-Adjusted Discontinuous
+            # Galerkin Methods: Curvilinear Meshes." Siam Journal on Scientific Computing
+            in_inv_op = op.RefInverseMassOperator(expr.op.dd_in, expr.op.dd_in)
+            trans_mass_op = op.RefMassOperator(expr.op.dd_in, expr.op.dd_out)
+            out_inv_op = op.RefInverseMassOperator(expr.op.dd_out, expr.op.dd_out)
+            return out_inv_op(trans_mass_op(in_inv_op(self.rec(expr.field))/jac_in))
 
         elif isinstance(expr.op, op.FaceMassOperator):
             jac_in_surf = sym.area_element(self.ambient_dim, self.dim - 1,
@@ -610,6 +616,9 @@ class StringifyMapper(pymbolic.mapper.stringifier.StringifyMapper):
     def map_nodal_min(self, expr, enclosing_prec):
         return "NodalMin" + self._format_op_dd(expr)
 
+    def map_elementwise_max(self, expr, enclosing_prec):
+        return "ElementwiseMax" + self._format_op_dd(expr)
+
     # }}}
 
     # {{{ global differentiation
@@ -718,7 +727,13 @@ class NoCSEStringifyMapper(StringifyMapper):
 # {{{ quadrature support
 
 class QuadratureCheckerAndRemover(CSECachingMapperMixin, IdentityMapper):
-    """Checks whether all quadratu
+    """ Checks that all quadrature tags exist in qttg (quad_tag_to_group_factory), 
+    and then removes all quadrature changes that are None in qttg.
+
+    For example: if you have an expression that interpolates to a QTAG called 
+    "product", and the discretization is made with qttg = {"product": None},
+    then all operators will be changed to ignore "product" QTAG.
+
     """
 
     def __init__(self, quad_tag_to_group_factory):
@@ -729,6 +744,7 @@ class QuadratureCheckerAndRemover(CSECachingMapperMixin, IdentityMapper):
     map_common_subexpression_uncached = \
             IdentityMapper.map_common_subexpression
 
+    # 
     def _process_dd(self, dd, location_descr):
         from grudge.symbolic.primitives import DOFDesc, QTAG_NONE
         if dd.quadrature_tag is not QTAG_NONE:
@@ -779,6 +795,58 @@ class QuadratureCheckerAndRemover(CSECachingMapperMixin, IdentityMapper):
         dd = self._process_dd(expr.dd, location_descr=type(expr).__name__)
         return type(expr)(expr.axis, dd)
 
+class QuadratureRemover(CSECachingMapperMixin, IdentityMapper):
+    """ Removes all quadrature use
+    """
+
+    def __init__(self):
+        IdentityMapper.__init__(self)
+        CSECachingMapperMixin.__init__(self)
+
+    map_common_subexpression_uncached = \
+            IdentityMapper.map_common_subexpression
+
+    def _process_dd(self, dd):
+        from grudge.symbolic.primitives import DOFDesc, QTAG_NONE
+        return DOFDesc(dd.domain_tag, QTAG_NONE)
+
+    def map_operator_binding(self, expr):
+        dd_in_orig = dd_in = expr.op.dd_in
+        dd_out_orig = dd_out = expr.op.dd_out
+
+        dd_in = self._process_dd(dd_in)
+        dd_out = self._process_dd(dd_out)
+
+        if dd_in_orig == dd_in and dd_out_orig == dd_out:
+            # unchanged
+            return IdentityMapper.map_operator_binding(self, expr)
+
+        import grudge.symbolic.operators as op
+        # changed
+
+        if dd_in == dd_out and isinstance(expr.op, op.InterpolationOperator):
+            # This used to be to-quad interpolation and has become a no-op.
+            # Remove it.
+            return self.rec(expr.field)
+
+        if isinstance(expr.op, op.StiffnessTOperator):
+            new_op = type(expr.op)(expr.op.xyz_axis, dd_in, dd_out)
+        elif isinstance(expr.op, (op.FaceMassOperator, op.InterpolationOperator)):
+            new_op = type(expr.op)(dd_in, dd_out)
+        else:
+            raise NotImplementedError("do not know how to modify dd_in and dd_out "
+                    "in %s" % type(expr.op).__name__)
+
+        return new_op(self.rec(expr.field))
+
+    #def map_ones(self, expr):
+        #dd = self._process_dd(expr.dd, location_descr=type(expr).__name__)
+        #return type(expr)(dd)
+
+    #def map_node_coordinate_component(self, expr):
+        #dd = self._process_dd(expr.dd, location_descr=type(expr).__name__)
+        #return type(expr)(expr.axis, dd)
+
 # }}}
 
 
diff --git a/grudge/symbolic/operators.py b/grudge/symbolic/operators.py
index 194b071614dba86a7ff150c075e5a3c987f23a4b..31dc316693ac02cece0e1bb8b893e4f985745ddd 100644
--- a/grudge/symbolic/operators.py
+++ b/grudge/symbolic/operators.py
@@ -57,6 +57,10 @@ class Operator(pymbolic.primitives.Expression):
     def __init__(self, dd_in, dd_out):
         self.dd_in = _sym().as_dofdesc(dd_in)
         self.dd_out = _sym().as_dofdesc(dd_out)
+        #TODO: remove this. According to Python docs this creates weird circular references
+        # And a bad time for everyone
+        #import inspect
+        #self.stack = inspect.stack()
 
     def stringifier(self):
         from grudge.symbolic.mappers import StringifyMapper
@@ -113,7 +117,7 @@ interp = InterpolationOperator
 
 class ElementwiseMaxOperator(Operator):
     mapper_method = intern("map_elementwise_max")
-
+    
 
 # {{{ nodal reduction: sum, integral, max
 
diff --git a/grudge/symbolic/primitives.py b/grudge/symbolic/primitives.py
index 00f7b8d386af6130810cf900bbaaebb0027f5adc..a29e481aa4e8253da7a7ae13e18234343507c2f3 100644
--- a/grudge/symbolic/primitives.py
+++ b/grudge/symbolic/primitives.py
@@ -198,6 +198,13 @@ class DOFDesc(object):
         self.domain_tag = domain_tag
         self.quadrature_tag = quadrature_tag
 
+
+        #TODO: remove this. According to Python docs this creates weird circular references
+        # And a bad time for everyone
+        #import inspect
+        #self.stack = inspect.stack()
+
+
     def is_scalar(self):
         return self.domain_tag is DTAG_SCALAR
 
@@ -311,6 +318,7 @@ class Variable(HasDOFDesc, ExpressionBase, pymbolic.primitives.Variable):
         HasDOFDesc.__init__(self, dd)
         pymbolic.primitives.Variable.__init__(self, name)
 
+
     def __getinitargs__(self):
         return (self.name, self.dd,)
 
diff --git a/test/test_grudge.py b/test/test_grudge.py
index 34532bfc80371e29bde28ea9674c48e36ab889aa..3976f246c6f1714a10730695df4cc52850657ed6 100644
--- a/test/test_grudge.py
+++ b/test/test_grudge.py
@@ -333,6 +333,80 @@ def test_convergence_advec(ctx_factory, mesh_name, mesh_pars, op_type, flux_type
     assert eoc_rec.order_estimate() > order
 
 
+@pytest.mark.parametrize("order", [3, 4, 5])
+def test_convergence_curvi_advection(ctx_factory,  order):
+    """Test whether a curvilinear mesh converges"""
+
+    cl_ctx = cl.create_some_context()
+    queue = cl.CommandQueue(cl_ctx)
+
+    from pytools.convergence import EOCRecorder
+    eoc_rec = EOCRecorder()
+
+    dims = 2
+    ns = [4, 6, 8]
+    for n in ns:
+
+        from meshmode.mesh.generation import generate_warped_rect_mesh
+        mesh = generate_warped_rect_mesh(dims, order, n)
+
+        discr = DGDiscretizationWithBoundaries(cl_ctx, mesh, order=order)
+
+        advec_v = np.array([0.1, 0.1])
+        norm_v = la.norm(advec_v)
+
+        def f(x):
+            return sym.sin(3*x)
+
+        def u_analytic(x):
+            return f(-np.dot(advec_v, x)/norm_v+sym.var("t", sym.DD_SCALAR)*norm_v)
+
+        from grudge.models.advection import WeakAdvectionOperator
+        analytic_sol = bind(discr, u_analytic(sym.nodes(dims)))
+
+        analytic_init = bind(discr, u_analytic(sym.nodes(dims)))
+
+        fields = analytic_init(queue, t=0)
+
+        op = WeakAdvectionOperator(advec_v,
+            inflow_u=u_analytic(sym.nodes(dims, sym.BTAG_ALL)),
+            flux_type="central")
+
+        bound_op = bind(discr, op.sym_operator())
+
+        def rhs(t, u):
+            return bound_op(queue, t=t, u=u)
+
+        final_t = 0.1
+
+        h = 1.0/n
+        dt_factor = 4
+        dt = dt_factor * h/order**2
+        nsteps = (final_t // dt) + 1
+        dt = final_t/nsteps + 1e-15
+
+        from grudge.shortcuts import set_up_rk4
+        dt_stepper = set_up_rk4("u", dt, fields, rhs)
+
+        norm = bind(discr, sym.norm(2, sym.var("u")))
+
+        step = 0
+        for event in dt_stepper.run(t_end=final_t):
+            if isinstance(event, dt_stepper.StateComputed):
+                assert event.component_id == "u"
+                esc = event.state_component
+                step += 1
+
+        sol = analytic_sol(queue, t=step * dt)
+        total_error = norm(queue, u=(esc - sol)) / norm(queue, u=sol)
+        eoc_rec.add_data_point(1.0/n, total_error)
+
+    print(eoc_rec.pretty_print(abscissa_label="h",
+            error_label="L2 Error"))
+
+    assert eoc_rec.order_estimate() > order
+
+
 @pytest.mark.parametrize("order", [3, 4, 5])
 def test_convergence_maxwell(ctx_factory,  order):
     """Test whether 3D maxwells actually converges"""