examples/wave/var2.py

+"""Minimal example of a grudge driver."""
+
+from __future__ import division, print_function
+
+__copyright__ = "Copyright (C) 2015 Andreas Kloeckner"
+
+__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
+from grudge.shortcuts import set_up_rk4
+from grudge import sym, bind, Discretization
+
+
+def main(write_output=True, order=4):
+    cl_ctx = cl.create_some_context()
+    queue = cl.CommandQueue(cl_ctx)
+
+    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=(8,)*dims)
+
+    discr = Discretization(cl_ctx, mesh, order=order)
+
+    source_center = np.array([0.1, 0.22, 0.33])[:mesh.dim]
+    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")
+    c =    sym.If(sym.Comparison(
+                np.dot(sym_x, sym_x), "<", 0.2),
+                -0.1, -0.2)
+    #c = np.dot(sym_x,sym_x)
+
+    from grudge.models.wave import VariableCoefficientWeakWaveOperator
+    from meshmode.mesh import BTAG_ALL, BTAG_NONE
+    op = VariableCoefficientWeakWaveOperator(c,
+            discr.dim,
+            source_f=(
+                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="central")
+
+    queue = cl.CommandQueue(discr.cl_context)
+    from pytools.obj_array import join_fields
+    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)
+
+    if mesh.dim == 2:
+        dt = 0.04
+    elif mesh.dim == 3:
+        dt = 0.02
+
+    dt_stepper = set_up_rk4("w", dt, fields, rhs)
+
+    final_t = 10
+    nsteps = int(final_t/dt)
+    print("dt=%g nsteps=%d" % (dt, nsteps))
+
+    from grudge.shortcuts import make_visualizer
+    vis = make_visualizer(discr, vis_order=order)
+
+    step = 0
+
+    norm = bind(discr, sym.norm(2, sym.var("u")))
+
+    from time import time
+    t_last_step = time()
+
+    for event in dt_stepper.run(t_end=final_t):
+        if isinstance(event, dt_stepper.StateComputed):
+            assert event.component_id == "w"
+
+            step += 1
+
+            print(step, event.t, norm(queue, u=event.state_component[0]),
+                    time()-t_last_step)
+            if step % 10 == 0:
+                vis.write_vtk_file("fld-%04d.vtu" % step,
+                        [
+                            ("u", event.state_component[0]),
+                            ("v", event.state_component[1:]),
+                            ])
+            t_last_step = time()
+
+
+if __name__ == "__main__":
+    main()

examples/wave/weak.py

+"""Minimal example of a grudge driver."""
+
+from __future__ import division, print_function
+
+__copyright__ = "Copyright (C) 2015 Andreas Kloeckner"
+
+__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
+from grudge.shortcuts import set_up_rk4
+from grudge import sym, bind, Discretization
+
+
+def main(write_output=True, order=4):
+    cl_ctx = cl.create_some_context()
+    queue = cl.CommandQueue(cl_ctx)
+
+    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=(8,)*dims)
+
+    discr = Discretization(cl_ctx, mesh, order=order)
+
+    source_center = np.array([0.1, 0.22, 0.33])[:mesh.dim]
+    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 grudge.models.wave import WeakWaveOperator
+    from meshmode.mesh import BTAG_ALL, BTAG_NONE
+    op = WeakWaveOperator(-0.1, discr.dim,
+            source_f=(
+                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")
+
+    queue = cl.CommandQueue(discr.cl_context)
+    from pytools.obj_array import join_fields
+    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)
+
+    if mesh.dim == 2:
+        dt = 0.04
+    elif mesh.dim == 3:
+        dt = 0.02
+
+    dt_stepper = set_up_rk4("w", dt, fields, rhs)
+
+    final_t = 10
+    nsteps = int(final_t/dt)
+    print("dt=%g nsteps=%d" % (dt, nsteps))
+
+    from grudge.shortcuts import make_visualizer
+    vis = make_visualizer(discr, vis_order=order)
+
+    step = 0
+
+    norm = bind(discr, sym.norm(2, sym.var("u")))
+
+    from time import time
+    t_last_step = time()
+
+    for event in dt_stepper.run(t_end=final_t):
+        if isinstance(event, dt_stepper.StateComputed):
+            assert event.component_id == "w"
+
+            step += 1
+
+            print(step, event.t, norm(queue, u=event.state_component[0]),
+                    time()-t_last_step)
+            if step % 10 == 0:
+                vis.write_vtk_file("fld-%04d.vtu" % step,
+                        [
+                            ("u", event.state_component[0]),
+                            ("v", event.state_component[1:]),
+                            ])
+            t_last_step = time()
+
+
+if __name__ == "__main__":
+    main()

grudge/execution.py

@@ -119,6 +119,17 @@ class ExecutionMapper(mappers.Evaluator,
         then = self.rec(expr.then)
         else_ = self.rec(expr.else_)
 
+        if (type(then) is not pyopencl.array.Array):
+            fill_val = then
+            then = cl.array.empty(self.queue,bool_crit.shape, dtype=np.float32) #what type should this be lol?
+            then.fill(fill_val) 
+
+        if (type(else_) is not pyopencl.array.Array):
+            fill_val = else_
+            else_ = cl.array.empty(self.queue,bool_crit.shape, dtype=np.float32)
+            else_.fill(fill_val)
+
+
         result = cl.array.empty_like(then, queue=self.queue,
                 allocator=self.bound_op.allocator)
         cl.array.if_positive(bool_crit, then, else_, out=result,

grudge/models/wave.py

@@ -180,187 +180,325 @@ class StrongWaveOperator(HyperbolicOperator):
     def max_eigenvalue(self, t, fields=None, discr=None):
         return abs(self.c)
 
-# }}}
-
-
-# {{{ variable-velocity
-
-class VariableVelocityStrongWaveOperator(HyperbolicOperator):
-    r"""This operator discretizes the wave equation
-    :math:`\partial_t^2 u = c^2 \Delta u`.
+class WeakWaveOperator(HyperbolicOperator):
+    """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 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.
 
-        \partial_t v - c \nabla u = 0
+    :math:`c` is assumed to be constant across all space.
     """
 
-    def __init__(
-            self, c, ambient_dim, source=0,
+    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,
-            time_sign=1,
-            diffusion_coeff=None,
-            diffusion_scheme=CentralSecondDerivative()
-            ):
-        """*c* and *source* are optemplate expressions.
-        """
+            radiation_tag=BTAG_NONE):
         assert isinstance(ambient_dim, int)
 
         self.c = c
-        self.time_sign = time_sign
         self.ambient_dim = ambient_dim
-        self.source = source
+        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.flux_type = flux_type
+        self.dirichlet_bc_f = dirichlet_bc_f
 
-        self.diffusion_coeff = diffusion_coeff
-        self.diffusion_scheme = diffusion_scheme
+        self.flux_type = flux_type
 
-    # {{{ flux ----------------------------------------------------------------
-    def flux(self, w, c_vol):
+    def flux(self, w):
         u = w[0]
         v = w[1:]
-        normal = sym.normal(w.tag, self.ambient_dim)
+        normal = sym.normal(w.dd, self.ambient_dim)
 
-        c = sym.RestrictToBoundary(w.tag)(c_vol)
+        flux_weak = join_fields(
+                np.dot(v.avg, normal),
+                u.avg * normal)
 
-        flux = self.time_sign*1/2*join_fields(
-                c.ext * np.dot(v.ext, normal)
-                - c.int * np.dot(v.int, normal),
-                normal*(c.ext*u.ext - c.int*u.int))
 
         if self.flux_type == "central":
-            pass
+            return -self.c*flux_weak
         elif self.flux_type == "upwind":
-            flux += join_fields(
-                    c.ext*u.ext - c.int*u.int,
-                    c.ext*normal*np.dot(normal, v.ext)
-                    - c.int*normal*np.dot(normal, v.int)
-                    )
+            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
+
+    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.interp("vol", self.dirichlet_tag)(u))
+        dir_v = sym.cse(sym.interp("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 = join_fields(2*dir_g - dir_u, dir_v)
         else:
-            raise ValueError("invalid flux type '%s'" % self.flux_type)
+            dir_bc = join_fields(-dir_u, dir_v)
+
+        dir_bc = sym.cse(dir_bc, "dir_bc")
+
+        # neumann BCs ---------------------------------------------------------
+        neu_u = sym.cse(sym.interp("vol", self.neumann_tag)(u))
+        neu_v = sym.cse(sym.interp("vol", self.neumann_tag)(v))
+        neu_bc = sym.cse(join_fields(neu_u, -neu_v), "neu_bc")
+
+        # radiation BCs -------------------------------------------------------
+        rad_normal = sym.normal(self.radiation_tag, d)
+
+        rad_u = sym.cse(sym.interp("vol", self.radiation_tag)(u))
+        rad_v = sym.cse(sym.interp("vol", self.radiation_tag)(v))
 
-        return flux
+        rad_bc = sym.cse(join_fields(
+                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 bind_characteristic_velocity(self, discr):
-        from grudge.symbolic.operators import ElementwiseMaxOperator
+        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()( 
+                join_fields(
+                    -self.c*np.dot(sym.stiffness_t(self.ambient_dim), v),
+                    -self.c*(sym.stiffness_t(self.ambient_dim)*u)
+                    )
 
-        compiled = discr.compile(ElementwiseMaxOperator()(self.c))
 
-        def do(t, w, **extra_context):
-            return compiled(t=t, w=w, **extra_context)
+                - 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))
 
-        return do
+                    ) )
 
-    def sym_operator(self, with_sensor=False):
-        d = self.ambient_dim
+        result[0] += self.source_f
 
-        w = sym.make_sym_array("w", d+1)
-        u = w[0]
-        v = w[1:]
+        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)
 
-        flux_w = join_fields(self.c, w)
 
-        # {{{ boundary conditions
+# }}}
 
-        # Dirichlet
-        dir_c = sym.RestrictToBoundary(self.dirichlet_tag) * self.c
-        dir_u = sym.RestrictToBoundary(self.dirichlet_tag) * u
-        dir_v = sym.RestrictToBoundary(self.dirichlet_tag) * v
 
-        dir_bc = join_fields(dir_c, -dir_u, dir_v)
+# {{{ variable-velocity
 
-        # Neumann
-        neu_c = sym.RestrictToBoundary(self.neumann_tag) * self.c
-        neu_u = sym.RestrictToBoundary(self.neumann_tag) * u
-        neu_v = sym.RestrictToBoundary(self.neumann_tag) * v
+class VariableCoefficientWeakWaveOperator(HyperbolicOperator):
+    """This operator discretizes the wave equation
+    :math:`\\partial_t^2 u = c^2 \\Delta u`.
 
-        neu_bc = join_fields(neu_c, neu_u, -neu_v)
+    To be precise, we discretize the hyperbolic system
 
-        # Radiation
-        rad_normal = sym.make_normal(self.radiation_tag, d)
+    .. math::
 
-        rad_c = sym.RestrictToBoundary(self.radiation_tag) * self.c
-        rad_u = sym.RestrictToBoundary(self.radiation_tag) * u
-        rad_v = sym.RestrictToBoundary(self.radiation_tag) * v
+        \partial_t u - c \\nabla \\cdot v = 0
 
-        rad_bc = join_fields(
-                rad_c,
-                0.5*(rad_u - self.time_sign*np.dot(rad_normal, rad_v)),
-                0.5*rad_normal*(np.dot(rad_normal, rad_v) - self.time_sign*rad_u)
-                )
+        \partial_t v - c \\nabla u = 0
 
-        # }}}
+    The sign of :math:`v` determines whether we discretize the forward or the
+    backward wave equation.
 
-        # {{{ diffusion -------------------------------------------------------
-        from pytools.obj_array import with_object_array_or_scalar
+    :math:`c` is assumed to be constant across all space.
+    """
 
-        def make_diffusion(arg):
-            if with_sensor or (
-                    self.diffusion_coeff is not None and self.diffusion_coeff != 0):
-                if self.diffusion_coeff is None:
-                    diffusion_coeff = 0
-                else:
-                    diffusion_coeff = self.diffusion_coeff
+    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)
 
-                if with_sensor:
-                    diffusion_coeff += sym.Field("sensor")
+        self.c = c
+        self.ambient_dim = ambient_dim
+        self.source_f = source_f
 
-                from grudge.second_order import SecondDerivativeTarget
+        self.sign = sym.If(sym.Comparison(
+                            self.c, ">", 0),
+                                            1, -1)
 
-                # strong_form here allows the reuse the value of grad u.
-                grad_tgt = SecondDerivativeTarget(
-                        self.ambient_dim, strong_form=True,
-                        operand=arg)
+        self.dirichlet_tag = dirichlet_tag
+        self.neumann_tag = neumann_tag
+        self.radiation_tag = radiation_tag
 
-                self.diffusion_scheme.grad(grad_tgt, bc_getter=None,
-                        dirichlet_tags=[], neumann_tags=[])
+        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)
+
+        surf_c = sym.interp("vol", u.dd)(self.c)
+
+        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
 
-                div_tgt = SecondDerivativeTarget(
-                        self.ambient_dim, strong_form=False,
-                        operand=diffusion_coeff*grad_tgt.minv_all)
+        #return self.c*flux_strong
 
-                self.diffusion_scheme.div(div_tgt,
-                        bc_getter=None,
-                        dirichlet_tags=[], neumann_tags=[])
+    def sym_operator(self):
+        d = self.ambient_dim
 
-                return div_tgt.minv_all
-            else:
-                return 0
+        w = sym.make_sym_array("w", d+1)
+        u = w[0]
+        v = w[1:]
 
-        # }}}
+        # boundary conditions -------------------------------------------------
+
+        # dirichlet BCs -------------------------------------------------------
+        dir_u = sym.cse(sym.interp("vol", self.dirichlet_tag)(u))
+        dir_v = sym.cse(sym.interp("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 = join_fields(2*dir_g - dir_u, dir_v)
+        else:
+            dir_bc = join_fields(-dir_u, dir_v)
+
+        dir_bc = sym.cse(dir_bc, "dir_bc")
+
+        # neumann BCs ---------------------------------------------------------
+        neu_u = sym.cse(sym.interp("vol", self.neumann_tag)(u))
+        neu_v = sym.cse(sym.interp("vol", self.neumann_tag)(v))
+        neu_bc = sym.cse(join_fields(neu_u, -neu_v), "neu_bc")
+
+        # radiation BCs -------------------------------------------------------
+        rad_normal = sym.normal(self.radiation_tag, d)
+
+        rad_u = sym.cse(sym.interp("vol", self.radiation_tag)(u))
+        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)
+                ), "rad_bc")
 
         # entire operator -----------------------------------------------------
         nabla = sym.nabla(d)
-        flux_op = sym.get_flux_operator(self.flux())
 
-        return (
-                - join_fields(
-                    - self.time_sign*self.c*np.dot(nabla, v) - make_diffusion(u)
-                    + self.source,
+        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))
 
-                    -self.time_sign*self.c*(nabla*u) - with_object_array_or_scalar(
-                        make_diffusion, v)
+                    #)# )
+
+        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.InverseMassOperator() * (
-                    flux_op(flux_w)
-                    + flux_op(sym.BoundaryPair(flux_w, dir_bc, self.dirichlet_tag))
-                    + flux_op(sym.BoundaryPair(flux_w, neu_bc, self.neumann_tag))
-                    + flux_op(sym.BoundaryPair(flux_w, rad_bc, self.radiation_tag))
-                    ))
+
+
+                - 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
@@ -369,9 +507,9 @@ class VariableVelocityStrongWaveOperator(HyperbolicOperator):
             self.neumann_tag,
             self.radiation_tag])
 
-    def max_eigenvalue_expr(self):
-        import grudge.symbolic as sym
-        return sym.NodalMax()(sym.CFunction("fabs")(self.c))
+    def max_eigenvalue(self, t, fields=None, discr=None):
+        return abs(self.c)
+
 
 # }}}