diff --git a/examples/wave/var-propagation-speed.py b/examples/wave/var-propagation-speed.py
index 24719059263217f760745f84afa2851668775450..cb5d3e9ed00b355f671fc810058bf4b9e2898e0c 100644
--- a/examples/wave/var-propagation-speed.py
+++ b/examples/wave/var-propagation-speed.py
@@ -51,10 +51,9 @@ def main(write_output=True, order=4):
     sym_x = sym.nodes(mesh.dim)
     sym_source_center_dist = sym_x - source_center
     sym_t = sym.ScalarVariable("t")
-    c =    sym.If(sym.Comparison(
-                np.dot(sym_x, sym_x), "<", 0.2),
-                -0.1, -0.2)
-    #c = np.dot(sym_x,sym_x)
+    c = sym.If(sym.Comparison(
+                np.dot(sym_x, sym_x), "<", 0.15),
+                np.float32(-0.1), np.float32(-0.2))
 
     from grudge.models.wave import VariableCoefficientWeakWaveOperator
     from meshmode.mesh import BTAG_ALL, BTAG_NONE
@@ -75,15 +74,9 @@ def main(write_output=True, order=4):
     fields = join_fields(discr.zeros(queue),
             [discr.zeros(queue) for i in range(discr.dim)])
 
-    # FIXME
-    #dt = op.estimate_rk4_timestep(discr, fields=fields)
-
     op.check_bc_coverage(mesh)
 
-    # print(sym.pretty(op.sym_operator()))
     bound_op = bind(discr, op.sym_operator())
-    print(bound_op)
-    # 1/0
 
     def rhs(t, w):
         return bound_op(queue, t=t, w=w)
diff --git a/grudge/execution.py b/grudge/execution.py
index f78c345f3c87f4c495afa3685b15680c06fae0cd..e3646b7361d8b1606a1a234fbe6089a51127f8b7 100644
--- a/grudge/execution.py
+++ b/grudge/execution.py
@@ -120,30 +120,32 @@ class ExecutionMapper(mappers.Evaluator,
         else_ = self.rec(expr.else_)
 
         import pymbolic.primitives as p
-
         var = p.Variable
 
-        i = var("i") #sym.Variable("i")
+        i = var("i")
         if isinstance(then,  pyopencl.array.Array):
-            sym_then = var("a")[i] #sym.Variable("a")[i]
+            sym_then = var("a")[i]
+        elif isinstance(then,  np.number):
+            sym_then = var("a")
         else:
-            sym_then = var("a") # sym.Variable("a")
-            then = np.float64(then)
+            raise TypeError(
+                "Expected parameter to be of type np.number or pyopencl.array.Array")
 
         if isinstance(else_,  pyopencl.array.Array):
-            sym_else = var("b")[i] # sym.Variable("b")[i]
+            sym_else = var("b")[i]
+        elif isinstance(then,  np.number):
+            sym_else = var("b")
         else:
-            sym_else = var("b") # sym.Variable("b")
-            else_ = np.float64(else_)
+            raise TypeError(
+                "Expected parameter to be of type np.number or pyopencl.array.Array")
 
         @memoize_in(self.bound_op, "map_if_knl")
         def knl():
             knl = lp.make_kernel(
-               "{[i]: 0<=i<n}",
+                "{[i]: 0<=i<n}",
                 [
-                    lp.Assignment(var("out")[i], p.If(var("crit")[i], sym_then, sym_else))
-
-                   #lp.Assignment(sym.Variable("out")[i], sym.If(sym.Variable("crit")[i], sym_then, sym_else))
+                    lp.Assignment(var("out")[i],
+                        p.If(var("crit")[i], sym_then, sym_else))
                 ])
             return lp.split_iname(knl, "i", 128, outer_tag="g.0", inner_tag="l.0")
 
diff --git a/grudge/models/wave.py b/grudge/models/wave.py
index 6ea1d923287761bba5efbb4998beda4a519e6318..019691257f57030a232ba0502ccc493624dc96d3 100644
--- a/grudge/models/wave.py
+++ b/grudge/models/wave.py
@@ -28,7 +28,6 @@ THE SOFTWARE.
 
 import numpy as np
 from grudge.models import HyperbolicOperator
-from grudge.models.second_order import CentralSecondDerivative
 from meshmode.mesh import BTAG_ALL, BTAG_NONE
 from grudge import sym
 from pytools.obj_array import join_fields
@@ -180,6 +179,7 @@ class StrongWaveOperator(HyperbolicOperator):
     def max_eigenvalue(self, t, fields=None, discr=None):
         return abs(self.c)
 
+
 class WeakWaveOperator(HyperbolicOperator):
     """This operator discretizes the wave equation
     :math:`\\partial_t^2 u = c^2 \\Delta u`.
@@ -232,26 +232,14 @@ class WeakWaveOperator(HyperbolicOperator):
                 np.dot(v.avg, normal),
                 u.avg * normal)
 
-
         if self.flux_type == "central":
             return -self.c*flux_weak
         elif self.flux_type == "upwind":
             return -self.c*(flux_weak + self.sign*join_fields(
                     0.5*(u.int-u.ext),
                     0.5*(normal * np.dot(normal, v.int-v.ext))))
-            # see doc/notes/grudge-notes.tm
-            # THis is terrible
-            #flux_weak -= self.sign*join_fields(
-                    #0.5*(u.int-u.ext),
-                    #0.5*(normal * np.dot(normal, v.int-v.ext)))
-        #else:
-            #raise ValueError("invalid flux type '%s'" % self.flux_type)
-
-        #flux_strong = join_fields(
-                #np.dot(v.int, normal),
-                #u.int * normal) - flux_weak
-
-        #return self.c*flux_strong
+        else:
+            raise ValueError("invalid flux type '%s'" % self.flux_type)
 
     def sym_operator(self):
         d = self.ambient_dim
@@ -295,39 +283,22 @@ class WeakWaveOperator(HyperbolicOperator):
                 ), "rad_bc")
 
         # entire operator -----------------------------------------------------
-        nabla = sym.nabla(d)
-
         def flux(pair):
             return sym.interp(pair.dd, "all_faces")(self.flux(pair))
 
-
-        #result = (
-                #- join_fields(
-                    #-self.c*np.dot(nabla, v),
-                    #-self.c*(nabla*u)
-                    #)
-                #+
-                #sym.InverseMassOperator()(
-                    #sym.FaceMassOperator()(
-                        #flux(sym.int_tpair(w))
-                        #+ flux(sym.bv_tpair(self.dirichlet_tag, w, dir_bc))
-                        #+ flux(sym.bv_tpair(self.neumann_tag, w, neu_bc))
-                        #+ flux(sym.bv_tpair(self.radiation_tag, w, rad_bc))
-                        #)))
-
-        result = sym.InverseMassOperator()( 
+        result = sym.InverseMassOperator()(
                 join_fields(
                     -self.c*np.dot(sym.stiffness_t(self.ambient_dim), v),
                     -self.c*(sym.stiffness_t(self.ambient_dim)*u)
                     )
 
 
-                - sym.FaceMassOperator()( flux(sym.int_tpair(w)) 
+                - sym.FaceMassOperator()(flux(sym.int_tpair(w))
                     + flux(sym.bv_tpair(self.dirichlet_tag, w, dir_bc))
                     + flux(sym.bv_tpair(self.neumann_tag, w, neu_bc))
                     + flux(sym.bv_tpair(self.radiation_tag, w, rad_bc))
 
-                    ) )
+                    ))
 
         result[0] += self.source_f
 
@@ -381,7 +352,7 @@ class VariableCoefficientWeakWaveOperator(HyperbolicOperator):
 
         self.sign = sym.If(sym.Comparison(
                             self.c, ">", 0),
-                                            1, -1)
+                                            np.int32(1), np.int32(-1))
 
         self.dirichlet_tag = dirichlet_tag
         self.neumann_tag = neumann_tag
@@ -401,24 +372,15 @@ class VariableCoefficientWeakWaveOperator(HyperbolicOperator):
         flux_weak = join_fields(
                 np.dot(v.avg, normal),
                 u.avg * normal)
-        return -surf_c*flux_weak
-
 
         if self.flux_type == "central":
             return -surf_c*flux_weak
-        #elif self.flux_type == "upwind":
-            # see doc/notes/grudge-notes.tm
-            #flux_weak -= self.sign*join_fields(
-                    #0.5*(u.int-u.ext),
-                    #0.5*(normal * np.dot(normal, v.int-v.ext)))
-        #else:
-            #raise ValueError("invalid flux type '%s'" % self.flux_type)
-
-        #flux_strong = join_fields(
-                #np.dot(v.int, normal),
-                #u.int * normal) - flux_weak
-
-        #return self.c*flux_strong
+        elif self.flux_type == "upwind":
+            return -self.c*(flux_weak + self.sign*join_fields(
+                    0.5*(u.int-u.ext),
+                    0.5*(normal * np.dot(normal, v.int-v.ext))))
+        else:
+            raise ValueError("invalid flux type '%s'" % self.flux_type)
 
     def sym_operator(self):
         d = self.ambient_dim
@@ -457,44 +419,29 @@ class VariableCoefficientWeakWaveOperator(HyperbolicOperator):
         rad_v = sym.cse(sym.interp("vol", self.radiation_tag)(v))
 
         rad_bc = sym.cse(join_fields(
-                0.5*(rad_u - sym.interp("vol",self.radiation_tag)(self.sign)*np.dot(rad_normal, rad_v)),
-                0.5*rad_normal*(np.dot(rad_normal, rad_v) -sym.interp("vol",self.radiation_tag)(self.sign)*rad_u)
+                0.5*(rad_u - sym.interp("vol", self.radiation_tag)(self.sign)
+                    * np.dot(rad_normal, rad_v)),
+                0.5*rad_normal*(np.dot(rad_normal, rad_v)
+                    - sym.interp("vol", self.radiation_tag)(self.sign)*rad_u)
                 ), "rad_bc")
 
         # entire operator -----------------------------------------------------
-        nabla = sym.nabla(d)
-
         def flux(pair):
             return sym.interp(pair.dd, "all_faces")(self.flux(pair))
 
-
-        #result = sym.InverseMassOperator()( 
-                #join_fields(
-                    #-self.c*np.dot(sym.stiffness_t(self.ambient_dim), v),
-                    #-self.c*(sym.stiffness_t(self.ambient_dim)*u)
-                    #)
-
-
-                #- sym.FaceMassOperator()( flux(sym.int_tpair(w)) 
-                    #+ flux(sym.bv_tpair(self.dirichlet_tag, w, dir_bc))
-                    #+ flux(sym.bv_tpair(self.neumann_tag, w, neu_bc))
-                    #+ flux(sym.bv_tpair(self.radiation_tag, w, rad_bc))
-
-                    #)# )
-
-        result = sym.InverseMassOperator()( 
+        result = sym.InverseMassOperator()(
                 join_fields(
                     -self.c*np.dot(sym.stiffness_t(self.ambient_dim), v),
                     -self.c*(sym.stiffness_t(self.ambient_dim)*u)
                     )
 
 
-                - sym.FaceMassOperator()( flux(sym.int_tpair(w)) 
+                - sym.FaceMassOperator()(flux(sym.int_tpair(w))
                     + flux(sym.bv_tpair(self.dirichlet_tag, w, dir_bc))
                     + flux(sym.bv_tpair(self.neumann_tag, w, neu_bc))
                     + flux(sym.bv_tpair(self.radiation_tag, w, rad_bc))
 
-                    ) )
+                    ))
 
         result[0] += self.source_f