diff --git a/examples/dagrt-fusion.py b/examples/dagrt-fusion.py
index a493059e642e24867c2b6658722cad837189d3ae..fd37f27f867bd9db0bd195371cd768a80a243f25 100755
--- a/examples/dagrt-fusion.py
+++ b/examples/dagrt-fusion.py
@@ -29,6 +29,25 @@ OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
 THE SOFTWARE.
 """
 
+# FIXME:
+# Results before https://github.com/inducer/grudge/pull/15 were better:
+#
+# Operator     | \parbox{1in}{\centering \% Memory Ops. Due to Scalar Assignments}
+# -------------+-------------------------------------------------------------------
+# 2D: Baseline | 51.1
+# 2D: Inlined  | 48.9
+# 3D: Baseline | 50.1
+# 3D: Inlined  | 48.6
+# INFO:__main__:Wrote '<stdout>'
+# ==== Scalar Assigment Inlining Impact ====
+# Operator     | Bytes Read | Bytes Written | Total      | \% of Baseline
+# -------------+------------+---------------+------------+----------------
+# 2D: Baseline | 9489600    | 3348000       | 12837600   | 100
+# 2D: Inlined  | 8949600    | 2808000       | 11757600   | 91.6
+# 3D: Baseline | 1745280000 | 505440000     | 2250720000 | 100
+# 3D: Inlined  | 1680480000 | 440640000     | 2121120000 | 94.2
+# INFO:__main__:Wrote '<stdout>'
+
 
 import contextlib
 import logging
@@ -413,66 +432,23 @@ class FusedRK4TimeStepper(RK4TimeStepperBase):
 
 # {{{ problem setup code
 
-def get_strong_wave_op_with_discr(actx, dims=2, order=4):
-    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)
-
-    logger.debug("%d elements", mesh.nelements)
-
-    discr = DGDiscretizationWithBoundaries(actx, mesh, order=order)
-
+def _get_source_term(dims):
     source_center = np.array([0.1, 0.22, 0.33])[:dims]
     source_width = 0.05
     source_omega = 3
 
-    sym_x = sym.nodes(mesh.dim)
+    sym_x = sym.nodes(dims)
     sym_source_center_dist = sym_x - source_center
     sym_t = sym.ScalarVariable("t")
 
-    from grudge.models.wave import StrongWaveOperator
-    from meshmode.mesh import BTAG_ALL, BTAG_NONE
-    op = StrongWaveOperator(-0.1, dims,
-            source_f=(
+    return (
                 sym.sin(source_omega*sym_t)
                 * sym.exp(
                     -np.dot(sym_source_center_dist, sym_source_center_dist)
-                    / source_width**2)),
-            dirichlet_tag=BTAG_NONE,
-            neumann_tag=BTAG_NONE,
-            radiation_tag=BTAG_ALL,
-            flux_type="upwind")
-
-    op.check_bc_coverage(mesh)
-
-    return (op.sym_operator(), discr)
-
-
-def dg_flux(c, tpair):
-    u = tpair[0]
-    v = tpair[1:]
-
-    dims = len(v)
-
-    normal = sym.normal(tpair.dd, dims)
-    flux_weak = flat_obj_array(
-            np.dot(v.avg, normal),
-            u.avg * normal)
-
-    flux_weak -= (1 if c > 0 else -1)*flat_obj_array(
-            0.5*(u.int-u.ext),
-            0.5*(normal * np.dot(normal, v.int-v.ext)))
-
-    flux_strong = flat_obj_array(
-            np.dot(v.int, normal),
-            u.int * normal) - flux_weak
-
-    return sym.project(tpair.dd, "all_faces")(c*flux_strong)
+                    / source_width**2))
 
 
-def get_strong_wave_op_with_discr_direct(actx, dims=2, order=4):
+def get_wave_op_with_discr(actx, dims=2, order=4):
     from meshmode.mesh.generation import generate_regular_rect_mesh
     mesh = generate_regular_rect_mesh(
             a=(-0.5,)*dims,
@@ -483,54 +459,21 @@ def get_strong_wave_op_with_discr_direct(actx, dims=2, order=4):
 
     discr = DGDiscretizationWithBoundaries(actx, mesh, order=order)
 
-    source_center = np.array([0.1, 0.22, 0.33])[:dims]
-    source_width = 0.05
-    source_omega = 3
-
-    sym_x = sym.nodes(mesh.dim)
-    sym_source_center_dist = sym_x - source_center
-    sym_t = sym.ScalarVariable("t")
-
-    from meshmode.mesh import BTAG_ALL
-
-    c = -0.1
-    sign = -1
-
-    w = sym.make_sym_array("w", dims+1)
-    u = w[0]
-    v = w[1:]
-
-    source_f = (
-                sym.sin(source_omega*sym_t)
-                * sym.exp(
-                    -np.dot(sym_source_center_dist, sym_source_center_dist)
-                    / source_width**2))
-
-    rad_normal = sym.normal(BTAG_ALL, dims)
-
-    rad_u = sym.cse(sym.project("vol", BTAG_ALL)(u))
-    rad_v = sym.cse(sym.project("vol", BTAG_ALL)(v))
-
-    rad_bc = sym.cse(flat_obj_array(
-            0.5*(rad_u - sign*np.dot(rad_normal, rad_v)),
-            0.5*rad_normal*(np.dot(rad_normal, rad_v) - sign*rad_u)
-            ), "rad_bc")
+    from grudge.models.wave import WeakWaveOperator
+    from meshmode.mesh import BTAG_ALL, BTAG_NONE
+    op = WeakWaveOperator(0.1, dims,
+            source_f=_get_source_term(dims),
+            dirichlet_tag=BTAG_NONE,
+            neumann_tag=BTAG_NONE,
+            radiation_tag=BTAG_ALL,
+            flux_type="upwind")
 
-    sym_operator = (
-            - flat_obj_array(
-                -c*np.dot(sym.nabla(dims), v) - source_f,
-                -c*(sym.nabla(dims)*u)
-                )
-            + sym.InverseMassOperator()(
-                sym.FaceMassOperator()(
-                    dg_flux(c, sym.int_tpair(w))
-                    + dg_flux(c, sym.bv_tpair(BTAG_ALL, w, rad_bc))
-                    )))
+    op.check_bc_coverage(mesh)
 
-    return (sym_operator, discr)
+    return (op.sym_operator(), discr)
 
 
-def get_strong_wave_component(state_component):
+def get_wave_component(state_component):
     return (state_component[0], state_component[1:])
 
 # }}}
@@ -545,10 +488,10 @@ def test_stepper_equivalence(ctx_factory, order=4):
 
     dims = 2
 
-    sym_operator, _ = get_strong_wave_op_with_discr(
-            actx, dims=dims, order=order)
-    sym_operator_direct, discr = get_strong_wave_op_with_discr_direct(
+    sym_operator, discr = get_wave_op_with_discr(
             actx, dims=dims, order=order)
+    #sym_operator_direct, discr = get_wave_op_with_discr_direct(
+    #        actx, dims=dims, order=order)
 
     if dims == 2:
         dt = 0.04
@@ -561,11 +504,11 @@ def test_stepper_equivalence(ctx_factory, order=4):
     bound_op = bind(discr, sym_operator)
 
     stepper = RK4TimeStepper(
-            discr, "w", bound_op, 1 + discr.dim, get_strong_wave_component)
+            discr, "w", bound_op, 1 + discr.dim, get_wave_component)
 
     fused_stepper = FusedRK4TimeStepper(
-            discr, "w", sym_operator_direct, 1 + discr.dim,
-            get_strong_wave_component)
+            discr, "w", sym_operator, 1 + discr.dim,
+            get_wave_component)
 
     t_start = 0
     t_end = 0.5
@@ -807,7 +750,7 @@ def test_assignment_memory_model(ctx_factory):
     queue = cl.CommandQueue(cl_ctx)
     actx = PyOpenCLArrayContext(queue)
 
-    _, discr = get_strong_wave_op_with_discr_direct(actx, dims=2, order=3)
+    _, discr = get_wave_op_with_discr(actx, dims=2, order=3)
 
     # Assignment instruction
     bound_op = bind(
@@ -837,7 +780,7 @@ def test_stepper_mem_ops(ctx_factory, use_fusion):
 
     dims = 2
 
-    sym_operator, discr = get_strong_wave_op_with_discr_direct(
+    sym_operator, discr = get_wave_op_with_discr(
             actx, dims=dims, order=3)
 
     t_start = 0
@@ -854,13 +797,13 @@ def test_stepper_mem_ops(ctx_factory, use_fusion):
 
         stepper = RK4TimeStepper(
                 discr, "w", bound_op, 1 + discr.dim,
-                get_strong_wave_component,
+                get_wave_component,
                 exec_mapper_factory=ExecutionMapperWithMemOpCounting)
 
     else:
         stepper = FusedRK4TimeStepper(
                 discr, "w", sym_operator, 1 + discr.dim,
-                get_strong_wave_component,
+                get_wave_component,
                 exec_mapper_factory=ExecutionMapperWithMemOpCounting)
 
     step = 0
@@ -1008,7 +951,7 @@ def test_stepper_timing(ctx_factory, use_fusion):
 
     dims = 3
 
-    sym_operator, discr = get_strong_wave_op_with_discr_direct(
+    sym_operator, discr = get_wave_op_with_discr(
             actx, dims=dims, order=3)
 
     t_start = 0
@@ -1025,13 +968,13 @@ def test_stepper_timing(ctx_factory, use_fusion):
 
         stepper = RK4TimeStepper(
                 discr, "w", bound_op, 1 + discr.dim,
-                get_strong_wave_component,
+                get_wave_component,
                 exec_mapper_factory=ExecutionMapperWithTiming)
 
     else:
         stepper = FusedRK4TimeStepper(
                 discr, "w", sym_operator, 1 + discr.dim,
-                get_strong_wave_component,
+                get_wave_component,
                 exec_mapper_factory=ExecutionMapperWithTiming)
 
     step = 0
@@ -1059,7 +1002,7 @@ def test_stepper_timing(ctx_factory, use_fusion):
 def get_example_stepper(actx, dims=2, order=3, use_fusion=True,
                         exec_mapper_factory=ExecutionMapper,
                         return_ic=False):
-    sym_operator, discr = get_strong_wave_op_with_discr_direct(
+    sym_operator, discr = get_wave_op_with_discr(
             actx, dims=dims, order=3)
 
     if not use_fusion:
@@ -1069,13 +1012,13 @@ def get_example_stepper(actx, dims=2, order=3, use_fusion=True,
 
         stepper = RK4TimeStepper(
                 discr, "w", bound_op, 1 + discr.dim,
-                get_strong_wave_component,
+                get_wave_component,
                 exec_mapper_factory=exec_mapper_factory)
 
     else:
         stepper = FusedRK4TimeStepper(
                 discr, "w", sym_operator, 1 + discr.dim,
-                get_strong_wave_component,
+                get_wave_component,
                 exec_mapper_factory=exec_mapper_factory)
 
     if return_ic:
@@ -1130,7 +1073,7 @@ def problem_stats(order=3):
     actx = PyOpenCLArrayContext(queue)
 
     with open_output_file("grudge-problem-stats.txt") as outf:
-        _, dg_discr_2d = get_strong_wave_op_with_discr_direct(
+        _, dg_discr_2d = get_wave_op_with_discr(
             actx, dims=2, order=order)
         print("Number of 2D elements:", dg_discr_2d.mesh.nelements, file=outf)
         vol_discr_2d = dg_discr_2d.discr_from_dd("vol")
@@ -1138,7 +1081,7 @@ def problem_stats(order=3):
         from pytools import one
         print("Number of DOFs per 2D element:", one(dofs_2d), file=outf)
 
-        _, dg_discr_3d = get_strong_wave_op_with_discr_direct(
+        _, dg_discr_3d = get_wave_op_with_discr(
             actx, dims=3, order=order)
         print("Number of 3D elements:", dg_discr_3d.mesh.nelements, file=outf)
         vol_discr_3d = dg_discr_3d.discr_from_dd("vol")
diff --git a/examples/wave/wave-min-mpi.py b/examples/wave/wave-min-mpi.py
index 95bcab7034cc1adc6a3feb30f700aa49d804ae55..3e0063d866942bcc921f583c9464ee6d78f33bd7 100644
--- a/examples/wave/wave-min-mpi.py
+++ b/examples/wave/wave-min-mpi.py
@@ -79,9 +79,9 @@ def main(write_output=True, order=4):
     sym_source_center_dist = sym_x - source_center
     sym_t = sym.ScalarVariable("t")
 
-    from grudge.models.wave import StrongWaveOperator
+    from grudge.models.wave import WeakWaveOperator
     from meshmode.mesh import BTAG_ALL, BTAG_NONE
-    op = StrongWaveOperator(-0.1, discr.dim,
+    op = WeakWaveOperator(0.1, discr.dim,
             source_f=(
                 sym.sin(source_omega*sym_t)
                 * sym.exp(
diff --git a/examples/wave/wave-min.py b/examples/wave/wave-min.py
index de6c9f92e75d449e6a42de80f7b6fa7270df8a15..e567251aa9631b1fe1da7e648208d91802c42ad7 100644
--- a/examples/wave/wave-min.py
+++ b/examples/wave/wave-min.py
@@ -61,9 +61,9 @@ def main(write_output=True, order=4):
     sym_source_center_dist = sym_x - source_center
     sym_t = sym.ScalarVariable("t")
 
-    from grudge.models.wave import StrongWaveOperator
+    from grudge.models.wave import WeakWaveOperator
     from meshmode.mesh import BTAG_ALL, BTAG_NONE
-    op = StrongWaveOperator(-0.1, discr.dim,
+    op = WeakWaveOperator(0.1, discr.dim,
             source_f=(
                 sym.sin(source_omega*sym_t)
                 * sym.exp(
diff --git a/grudge/models/wave.py b/grudge/models/wave.py
index 1be9abe2a7339870941fcbc8e973a7eea7137c3e..863b0112052f71853fd746a80e29385b5ebf187f 100644
--- a/grudge/models/wave.py
+++ b/grudge/models/wave.py
@@ -35,153 +35,6 @@ from pytools.obj_array import flat_obj_array
 
 # {{{ constant-velocity
 
-class StrongWaveOperator(HyperbolicOperator):
-    r"""This operator discretizes the wave equation
-    :math:`\partial_t^2 u = c^2 \Delta u`.
-
-    To be precise, we discretize the hyperbolic system
-
-    .. math::
-
-        \partial_t u - c \\nabla \\cdot v = 0
-
-        \partial_t v - c \\nabla u = 0
-
-    The sign of :math:`v` determines whether we discretize the forward or the
-    backward wave equation.
-
-    :math:`c` is assumed to be constant across all space.
-    """
-
-    def __init__(self, c, ambient_dim, source_f=0,
-            flux_type="upwind",
-            dirichlet_tag=BTAG_ALL,
-            dirichlet_bc_f=0,
-            neumann_tag=BTAG_NONE,
-            radiation_tag=BTAG_NONE):
-        assert isinstance(ambient_dim, int)
-
-        self.c = c
-        self.ambient_dim = ambient_dim
-        self.source_f = source_f
-
-        if self.c > 0:
-            self.sign = 1
-        else:
-            self.sign = -1
-
-        self.dirichlet_tag = dirichlet_tag
-        self.neumann_tag = neumann_tag
-        self.radiation_tag = radiation_tag
-
-        self.dirichlet_bc_f = dirichlet_bc_f
-
-        self.flux_type = flux_type
-
-    def flux(self, w):
-        u = w[0]
-        v = w[1:]
-        normal = sym.normal(w.dd, self.ambient_dim)
-
-        flux_weak = flat_obj_array(
-                np.dot(v.avg, normal),
-                u.avg * normal)
-
-        if self.flux_type == "central":
-            pass
-        elif self.flux_type == "upwind":
-            # see doc/notes/grudge-notes.tm
-            flux_weak -= self.sign*flat_obj_array(
-                    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 = flat_obj_array(
-                np.dot(v.int, normal),
-                u.int * normal) - flux_weak
-
-        return self.c*flux_strong
-
-    def sym_operator(self):
-        d = self.ambient_dim
-
-        w = sym.make_sym_array("w", d+1)
-        u = w[0]
-        v = w[1:]
-
-        # boundary conditions -------------------------------------------------
-
-        # dirichlet BCs -------------------------------------------------------
-        dir_u = sym.cse(sym.project("vol", self.dirichlet_tag)(u))
-        dir_v = sym.cse(sym.project("vol", self.dirichlet_tag)(v))
-        if self.dirichlet_bc_f:
-            # FIXME
-            from warnings import warn
-            warn("Inhomogeneous Dirichlet conditions on the wave equation "
-                    "are still having issues.")
-
-            dir_g = sym.Field("dir_bc_u")
-            dir_bc = flat_obj_array(2*dir_g - dir_u, dir_v)
-        else:
-            dir_bc = flat_obj_array(-dir_u, dir_v)
-
-        dir_bc = sym.cse(dir_bc, "dir_bc")
-
-        # neumann BCs ---------------------------------------------------------
-        neu_u = sym.cse(sym.project("vol", self.neumann_tag)(u))
-        neu_v = sym.cse(sym.project("vol", self.neumann_tag)(v))
-        neu_bc = sym.cse(flat_obj_array(neu_u, -neu_v), "neu_bc")
-
-        # radiation BCs -------------------------------------------------------
-        rad_normal = sym.normal(self.radiation_tag, d)
-
-        rad_u = sym.cse(sym.project("vol", self.radiation_tag)(u))
-        rad_v = sym.cse(sym.project("vol", self.radiation_tag)(v))
-
-        rad_bc = sym.cse(flat_obj_array(
-                0.5*(rad_u - self.sign*np.dot(rad_normal, rad_v)),
-                0.5*rad_normal*(np.dot(rad_normal, rad_v) - self.sign*rad_u)
-                ), "rad_bc")
-
-        # entire operator -----------------------------------------------------
-        nabla = sym.nabla(d)
-
-        def flux(pair):
-            return sym.project(pair.dd, "all_faces")(
-                    self.flux(pair))
-
-        result = (
-                - flat_obj_array(
-                    -self.c*np.dot(nabla, v),
-                    -self.c*(nabla*u)
-                    )
-                +  # noqa: W504
-                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[0] += self.source_f
-
-        return result
-
-    def check_bc_coverage(self, mesh):
-        from meshmode.mesh import check_bc_coverage
-        check_bc_coverage(mesh, [
-            self.dirichlet_tag,
-            self.neumann_tag,
-            self.radiation_tag])
-
-    def max_eigenvalue(self, t, fields=None, discr=None):
-        return abs(self.c)
-
-
 class WeakWaveOperator(HyperbolicOperator):
     r"""This operator discretizes the wave equation
     :math:`\partial_t^2 u = c^2 \Delta u`.
@@ -230,16 +83,16 @@ class WeakWaveOperator(HyperbolicOperator):
         v = w[1:]
         normal = sym.normal(w.dd, self.ambient_dim)
 
-        flux_weak = flat_obj_array(
+        central_flux_weak = -self.c*flat_obj_array(
                 np.dot(v.avg, normal),
                 u.avg * normal)
 
         if self.flux_type == "central":
-            return -self.c*flux_weak
+            return central_flux_weak
         elif self.flux_type == "upwind":
-            return -self.c*(flux_weak + self.sign*flat_obj_array(
-                    0.5*(u.int-u.ext),
-                    0.5*(normal * np.dot(normal, v.int-v.ext))))
+            return central_flux_weak - self.c*self.sign*flat_obj_array(
+                    0.5*(u.ext-u.int),
+                    0.5*(normal * np.dot(normal, v.ext-v.int)))
         else:
             raise ValueError("invalid flux type '%s'" % self.flux_type)
 
@@ -294,7 +147,6 @@ class WeakWaveOperator(HyperbolicOperator):
                     -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))
diff --git a/test/test_mpi_communication.py b/test/test_mpi_communication.py
index a5109737106279f0a76d147a9d77a4e35d2d6d84..91656bda4d1e123cc6a69eafe0c2f1a9b4eb1fa4 100644
--- a/test/test_mpi_communication.py
+++ b/test/test_mpi_communication.py
@@ -140,9 +140,9 @@ def mpi_communication_entrypoint():
     sym_source_center_dist = sym_x - source_center
     sym_t = sym.ScalarVariable("t")
 
-    from grudge.models.wave import StrongWaveOperator
+    from grudge.models.wave import WeakWaveOperator
     from meshmode.mesh import BTAG_ALL, BTAG_NONE
-    op = StrongWaveOperator(-0.1, vol_discr.dim,
+    op = WeakWaveOperator(0.1, vol_discr.dim,
             source_f=(
                 sym.sin(source_omega*sym_t)
                 * sym.exp(