examples/maxwell/analytic_solutions.py

@@ -20,9 +20,9 @@
 from __future__ import division
 from __future__ import absolute_import
 from __future__ import print_function
-from grudge.tools import \
-        cyl_bessel_j, \
-        cyl_bessel_j_prime
+#from grudge.tools import \
+        #cyl_bessel_j, \
+        #cyl_bessel_j_prime
 from math import sqrt, pi, sin, cos, atan2, acos
 import numpy
 import numpy.linalg as la

examples/maxwell/cavities.py

@@ -20,14 +20,17 @@ from __future__ import absolute_import
 from __future__ import print_function
 import numpy as np
 
-import logging
-logger = logging.getLogger(__name__)
+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, allow_features=None, flux_type_arg=1,
         bdry_flux_type_arg=None, extra_discr_args={}):
-    from grudge.mesh.generator import make_cylinder_mesh, make_box_mesh
-    from grudge.tools import EOCRecorder, to_obj_array
+    cl_ctx = cl.create_some_context()
+    queue = cl.CommandQueue(cl_ctx)
+
     from math import sqrt, pi  # noqa
     from analytic_solutions import (  # noqa
             RealPartAdapter,
@@ -36,21 +39,16 @@ def main(write_output=True, allow_features=None, flux_type_arg=1,
             CylindricalCavityMode,
             RectangularWaveguideMode,
             RectangularCavityMode)
-    from grudge.models.em import MaxwellOperator
 
-    logging.basicConfig(level=logging.DEBUG)
 
-    from grudge.backends import guess_run_context
-    rcon = guess_run_context(allow_features)
 
     epsilon0 = 8.8541878176e-12  # C**2 / (N m**2)
     mu0 = 4*pi*1e-7  # N/A**2.
     epsilon = 1*epsilon0
     mu = 1*mu0
 
-    eoc_rec = EOCRecorder()
 
-    cylindrical = False
+    cylindrical = False #True is unsupported because meshmode
     periodic = False
 
     if cylindrical:
@@ -62,36 +60,29 @@ def main(write_output=True, allow_features=None, flux_type_arg=1,
         # r_sol = CylindricalFieldAdapter(RealPartAdapter(mode))
         # c_sol = SplitComplexAdapter(CylindricalFieldAdapter(mode))
 
-        if rcon.is_head_rank:
-            mesh = make_cylinder_mesh(radius=R, height=d, max_volume=0.01)
+        #mesh = make_cylinder_mesh(radius=R, height=d, max_volume=0.01)
     else:
         if periodic:
-            mode = RectangularWaveguideMode(epsilon, mu, (3, 2, 1))
+            mode = RectangularWaveguideMode(epsilon, mu, (3, 2, 1)) #What is the point of this if it is overwritten?
             periodicity = (False, False, True)
         else:
             periodicity = None
         mode = RectangularCavityMode(epsilon, mu, (1, 2, 2))
 
-        if rcon.is_head_rank:
-            mesh = make_box_mesh(max_volume=0.001, periodicity=periodicity)
+        mesh = make_box_mesh(max_volume=0.001, periodicity=periodicity)
 
-    if rcon.is_head_rank:
-        mesh_data = rcon.distribute_mesh(mesh)
-    else:
-        mesh_data = rcon.receive_mesh()
 
     for order in [4, 5, 6]:
+        discr = Discretization(mesh_data, order=order,
+                tune_for=op.sym_operator(),
+                **extra_discr_args)
     #for order in [1,2,3,4,5,6]:
-        extra_discr_args.setdefault("debug", []).extend([
-            "cuda_no_plan", "cuda_dump_kernels"])
+        from grudge.models.em import MaxwellOperator
 
         op = MaxwellOperator(epsilon, mu,
                 flux_type=flux_type_arg,
                 bdry_flux_type=bdry_flux_type_arg)
 
-        discr = rcon.make_discretization(mesh_data, order=order,
-                tune_for=op.sym_operator(),
-                **extra_discr_args)
 
         from grudge.visualization import VtkVisualizer
         if write_output:
@@ -188,20 +179,6 @@ def main(write_output=True, allow_features=None, flux_type_arg=1,
             eoc_rec.add_data_point(order,
                     discr.norm(fields-get_true_field()))
 
-        finally:
-            if write_output:
-                vis.close()
-
-            logmgr.close()
-            discr.close()
-
-        if rcon.is_head_rank:
-            print()
-            print(eoc_rec.pretty_print("P.Deg.", "L2 Error"))
-
-        # }}}
-
-    assert eoc_rec.estimate_order_of_convergence()[0, 1] > 6
 
 
 # {{{ entry points for py.test

examples/maxwell/cavities2.py

+# grudge - the Hybrid'n'Easy DG Environment
+# 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/>.
+
+
+from __future__ import division
+from __future__ import absolute_import
+from __future__ import print_function
+import numpy as np
+
+import logging
+logger = logging.getLogger(__name__)
+
+
+def main(write_output=True, allow_features=None, flux_type_arg=1,
+        bdry_flux_type_arg=None, extra_discr_args={}):
+    from grudge.mesh.generator import make_cylinder_mesh, make_box_mesh
+    from grudge.tools import EOCRecorder, to_obj_array
+    from math import sqrt, pi  # noqa
+    from analytic_solutions import (  # noqa
+            RealPartAdapter,
+            SplitComplexAdapter,
+            CylindricalFieldAdapter,
+            CylindricalCavityMode,
+            RectangularWaveguideMode,
+            RectangularCavityMode)
+    from grudge.models.em import MaxwellOperator
+
+    logging.basicConfig(level=logging.DEBUG)
+
+    from grudge.backends import guess_run_context
+    rcon = guess_run_context(allow_features)
+
+    epsilon0 = 8.8541878176e-12  # C**2 / (N m**2)
+    mu0 = 4*pi*1e-7  # N/A**2.
+    epsilon = 1*epsilon0
+    mu = 1*mu0
+
+    eoc_rec = EOCRecorder()
+
+    cylindrical = False
+    periodic = False
+
+    if cylindrical:
+        R = 1
+        d = 2
+        mode = CylindricalCavityMode(m=1, n=1, p=1,
+                radius=R, height=d,
+                epsilon=epsilon, mu=mu)
+        # r_sol = CylindricalFieldAdapter(RealPartAdapter(mode))
+        # c_sol = SplitComplexAdapter(CylindricalFieldAdapter(mode))
+
+        if rcon.is_head_rank:
+            mesh = make_cylinder_mesh(radius=R, height=d, max_volume=0.01)
+    else:
+        if periodic:
+            mode = RectangularWaveguideMode(epsilon, mu, (3, 2, 1))
+            periodicity = (False, False, True)
+        else:
+            periodicity = None
+        mode = RectangularCavityMode(epsilon, mu, (1, 2, 2))
+
+        if rcon.is_head_rank:
+            mesh = make_box_mesh(max_volume=0.001, periodicity=periodicity)
+
+    if rcon.is_head_rank:
+        mesh_data = rcon.distribute_mesh(mesh)
+    else:
+        mesh_data = rcon.receive_mesh()
+
+    for order in [4, 5, 6]:
+    #for order in [1,2,3,4,5,6]:
+        extra_discr_args.setdefault("debug", []).extend([
+            "cuda_no_plan", "cuda_dump_kernels"])
+
+        op = MaxwellOperator(epsilon, mu,
+                flux_type=flux_type_arg,
+                bdry_flux_type=bdry_flux_type_arg)
+
+        discr = rcon.make_discretization(mesh_data, order=order,
+                tune_for=op.sym_operator(),
+                **extra_discr_args)
+
+        from grudge.visualization import VtkVisualizer
+        if write_output:
+            vis = VtkVisualizer(discr, rcon, "em-%d" % order)
+
+        mode.set_time(0)
+
+        def get_true_field():
+            return discr.convert_volume(
+                to_obj_array(mode(discr)
+                    .real.astype(discr.default_scalar_type).copy()),
+                kind=discr.compute_kind)
+
+        fields = get_true_field()
+
+        if rcon.is_head_rank:
+            print("---------------------------------------------")
+            print("order %d" % order)
+            print("---------------------------------------------")
+            print("#elements=", len(mesh.elements))
+
+        from grudge.timestep.runge_kutta import LSRK4TimeStepper
+        stepper = LSRK4TimeStepper(dtype=discr.default_scalar_type, rcon=rcon)
+        #from grudge.timestep.dumka3 import Dumka3TimeStepper
+        #stepper = Dumka3TimeStepper(3, dtype=discr.default_scalar_type, rcon=rcon)
+
+        # {{{ diagnostics setup
+
+        from pytools.log import LogManager, add_general_quantities, \
+                add_simulation_quantities, add_run_info
+
+        if write_output:
+            log_file_name = "maxwell-%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)
+
+        from pytools.log import IntervalTimer
+        vis_timer = IntervalTimer("t_vis", "Time spent visualizing")
+        logmgr.add_quantity(vis_timer)
+
+        from grudge.log import EMFieldGetter, add_em_quantities
+        field_getter = EMFieldGetter(discr, op, lambda: fields)
+        add_em_quantities(logmgr, op, field_getter)
+
+        logmgr.add_watches(
+                ["step.max", "t_sim.max",
+                    ("W_field", "W_el+W_mag"),
+                    "t_step.max"]
+                )
+
+        # }}}
+
+        # {{{ timestep loop
+
+        rhs = op.bind(discr)
+        final_time = 0.5e-9
+
+        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, fields=fields))
+
+            for step, t, dt in step_it:
+                if step % 50 == 0 and write_output:
+                    sub_timer = vis_timer.start_sub_timer()
+                    e, h = op.split_eh(fields)
+                    visf = vis.make_file("em-%d-%04d" % (order, step))
+                    vis.add_data(
+                            visf,
+                            [
+                                ("e",
+                                    discr.convert_volume(e, kind="numpy")),
+                                ("h",
+                                    discr.convert_volume(h, kind="numpy")),
+                                ],
+                            time=t, step=step)
+                    visf.close()
+                    sub_timer.stop().submit()
+
+                fields = stepper(fields, t, dt, rhs)
+
+            mode.set_time(final_time)
+
+            eoc_rec.add_data_point(order,
+                    discr.norm(fields-get_true_field()))
+
+        finally:
+            if write_output:
+                vis.close()
+
+            logmgr.close()
+            discr.close()
+
+        if rcon.is_head_rank:
+            print()
+            print(eoc_rec.pretty_print("P.Deg.", "L2 Error"))
+
+        # }}}
+
+    assert eoc_rec.estimate_order_of_convergence()[0, 1] > 6
+
+
+# {{{ entry points for py.test
+
+from pytools.test import mark_test
+
+
+@mark_test.long
+def test_maxwell_cavities():
+    main(write_output=False)
+
+
+@mark_test.long
+def test_maxwell_cavities_lf():
+    main(write_output=False, flux_type_arg="lf", bdry_flux_type_arg=1)
+
+
+@mark_test.mpi
+@mark_test.long
+def test_maxwell_cavities_mpi():
+    from pytools.mpi import run_with_mpi_ranks
+    run_with_mpi_ranks(__file__, 2, main,
+            write_output=False, allow_features=["mpi"])
+
+
+def test_cuda():
+    try:
+        from pycuda.tools import mark_cuda_test
+    except ImportError:
+        pass
+
+    yield "SP CUDA Maxwell", mark_cuda_test(
+            do_test_maxwell_cavities_cuda), np.float32
+    yield "DP CUDA Maxwell", mark_cuda_test(
+            do_test_maxwell_cavities_cuda), np.float64
+
+
+def do_test_maxwell_cavities_cuda(dtype):
+    import pytest  # noqa
+
+    main(write_output=False, allow_features=["cuda"],
+            extra_discr_args=dict(
+                debug=["cuda_no_plan"],
+                default_scalar_type=dtype,
+                ))
+
+# }}}
+
+
+# entry point -----------------------------------------------------------------
+
+if __name__ == "__main__":
+    from pytools.mpi import in_mpi_relaunch
+    if in_mpi_relaunch():
+        test_maxwell_cavities_mpi()
+    else:
+        do_test_maxwell_cavities_cuda(np.float32)

examples/maxwell/cavities3.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
+
+from analytic_solutions import (  # noqa
+        RectangularWaveguideMode,
+        RectangularCavityMode)
+
+def main(write_output=True, order=4):
+    cl_ctx = cl.create_some_context()
+    queue = cl.CommandQueue(cl_ctx)
+
+    dims = 3
+    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)
+
+
+    epsilon0 = 8.8541878176e-12  # C**2 / (N m**2)
+    mu0 = 4*np.pi*1e-7  # N/A**2.
+    epsilon = 1*epsilon0
+    mu = 1*mu0
+
+
+    from grudge.models.em import MaxwellOperator
+    from meshmode.mesh import BTAG_ALL, BTAG_NONE
+    op = MaxwellOperator(epsilon,mu,
+            flux_type="lf")
+    
+
+    queue = cl.CommandQueue(discr.cl_context)
+    from pytools.obj_array import join_fields
+    fields = RectangularCavityMode(epsilon, mu, (1, 2, 2))(discr)
+    #fields = join_fields(discr.zeros(queue),
+            #[discr.zeros(queue) for i in range(discr.dim*2 - 1)])
+
+    # 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,
+                        [
+                            ("e", event.state_component[0]),
+                            ("h", event.state_component[1:]),
+                            ])
+            t_last_step = time()
+
+
+if __name__ == "__main__":
+    main()

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="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()

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): # is instance, poarenthesis
+            fill_val = then
+            then = cl.array.empty(self.queue,bool_crit.shape, dtype=np.float64) #what type should this be?
+            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.float64)
+            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/em.py

@@ -30,8 +30,8 @@ THE SOFTWARE.
 from pytools import memoize_method
 
 from grudge.models import HyperbolicOperator
-from grudge.symbolic.primitives import make_common_subexpression as cse
 from meshmode.mesh import BTAG_ALL, BTAG_NONE
+from grudge import sym
 from pytools.obj_array import join_fields, make_obj_array
 
 # TODO: Check PML
@@ -106,17 +106,8 @@ class MaxwellOperator(HyperbolicOperator):
         self.current = current
         self.incident_bc_data = incident_bc
 
-    @property
-    def c(self):
-        from warnings import warn
-        warn("MaxwellOperator.c is deprecated", DeprecationWarning)
-        if not self.fixed_material:
-            raise RuntimeError("Cannot compute speed of light "
-                    "for non-constant material")
 
-        return 1/(self.mu*self.epsilon)**0.5
-
-    def flux(self, flux_type):
+    def flux(self, w):
         """The template for the numerical flux for variable coefficients.
 
         :param flux_type: can be in [0,1] for anything between central and upwind,
@@ -124,38 +115,23 @@ class MaxwellOperator(HyperbolicOperator):
 
         As per Hesthaven and Warburton page 433.
         """
-        from grudge.flux import (make_normal, FluxVectorPlaceholder,
-                FluxConstantPlaceholder)
 
-        normal = make_normal(self.dimensions)
+        normal = sym.normal(w.dd, self.dimensions)
 
         if self.fixed_material:
-            from grudge.tools import count_subset
-            w = FluxVectorPlaceholder(count_subset(self.get_eh_subset()))
-
             e, h = self.split_eh(w)
-            epsilon = FluxConstantPlaceholder(self.epsilon)
-            mu = FluxConstantPlaceholder(self.mu)
-
-        else:
-            from grudge.tools import count_subset
-            w = FluxVectorPlaceholder(count_subset(self.get_eh_subset())+2)
+            epsilon = self.epsilon
+            mu = self.mu
 
-            epsilon, mu, e, h = self.split_eps_mu_eh(w)
 
-        Z_int = (mu.int/epsilon.int)**0.5
+        Z_int = (mu/epsilon)**0.5
         Y_int = 1/Z_int
-        Z_ext = (mu.ext/epsilon.ext)**0.5
+        Z_ext = (mu/epsilon)**0.5
         Y_ext = 1/Z_ext
 
-        if flux_type == "lf":
+        if self.flux_type == "lf":
             if self.fixed_material:
                 max_c = (self.epsilon*self.mu)**(-0.5)
-            else:
-                from grudge.flux import Max
-                c_int = (epsilon.int*mu.int)**(-0.5)
-                c_ext = (epsilon.ext*mu.ext)**(-0.5)
-                max_c = Max(c_int, c_ext)  # noqa
 
             return join_fields(
                     # flux e,
@@ -170,42 +146,41 @@ class MaxwellOperator(HyperbolicOperator):
                         # multiplication by mu undoes material divisor below
                         #-max_c*(mu.int*h.int - mu.ext*h.ext)
                     ))
-        elif isinstance(flux_type, (int, float)):
+        elif isinstance(self.flux_type, (int, float)):
             # see doc/maxima/maxwell.mac
             return join_fields(
                     # flux e,
                     (
                         -1/(Z_int+Z_ext)*self.space_cross_h(normal,
                             Z_ext*(h.int-h.ext)
-                            - flux_type*self.space_cross_e(normal, e.int-e.ext))
+                            - self.flux_type*self.space_cross_e(normal, e.int-e.ext))
                         ),
                     # flux h
                     (
                         1/(Y_int + Y_ext)*self.space_cross_e(normal,
                             Y_ext*(e.int-e.ext)
-                            + flux_type*self.space_cross_h(normal, h.int-h.ext))
+                            + self.flux_type*self.space_cross_h(normal, h.int-h.ext))
                         ),
                     )
         else:
             raise ValueError("maxwell: invalid flux_type (%s)"
                     % self.flux_type)
 
-    def local_derivatives(self, w=None):
+    def local_derivatives(self, w):
         """Template for the spatial derivatives of the relevant components of
         :math:`E` and :math:`H`
         """
 
-        e, h = self.split_eh(self.field_placeholder(w))
+        e, h = self.split_eh(w)
 
+        nabla = sym.nabla(self.dimensions)
         def e_curl(field):
             return self.space_cross_e(nabla, field)
 
         def h_curl(field):
             return self.space_cross_h(nabla, field)
 
-        from grudge.symbolic import make_nabla
 
-        nabla = make_nabla(self.dimensions)
 
         # in conservation form: u_t + A u_x = 0
         return join_fields(
@@ -213,62 +188,44 @@ class MaxwellOperator(HyperbolicOperator):
                 e_curl(e)
                 )
 
-    def field_placeholder(self, w=None):
-        "A placeholder for E and H."
-        from grudge.tools import count_subset
-        fld_cnt = count_subset(self.get_eh_subset())
-        if w is None:
-            from grudge.symbolic import make_sym_vector
-            w = make_sym_vector("w", fld_cnt)
-
-        return w
-
-    def pec_bc(self, w=None):
+    def pec_bc(self, w):
         "Construct part of the flux operator template for PEC boundary conditions"
-        e, h = self.split_eh(self.field_placeholder(w))
+        e, h = self.split_eh(w)
+
+        pec_e = sym.cse(sym.interp("vol", self.pec_tag)(e))
+        pec_h = sym.cse(sym.interp("vol", self.pec_tag)(h))
 
-        from grudge.symbolic import RestrictToBoundary
-        pec_e = RestrictToBoundary(self.pec_tag)(e)
-        pec_h = RestrictToBoundary(self.pec_tag)(h)
 
         return join_fields(-pec_e, pec_h)
 
-    def pmc_bc(self, w=None):
+    def pmc_bc(self, w):
         "Construct part of the flux operator template for PMC boundary conditions"
-        e, h = self.split_eh(self.field_placeholder(w))
+        e, h = self.split_eh(w)
 
-        from grudge.symbolic import RestrictToBoundary
-        pmc_e = RestrictToBoundary(self.pmc_tag)(e)
-        pmc_h = RestrictToBoundary(self.pmc_tag)(h)
+        pmc_e = sym.cse(sym.interp("vol", self.pmc_tag)(e))
+        pmc_h = sym.cse(sym.interp("vol", self.pmc_tag)(h))
 
         return join_fields(pmc_e, -pmc_h)
 
-    def absorbing_bc(self, w=None):
+    def absorbing_bc(self, w):
         """Construct part of the flux operator template for 1st order
         absorbing boundary conditions.
         """
 
-        from grudge.symbolic import normal
-        absorb_normal = normal(self.absorb_tag, self.dimensions)
+        absorb_normal = sym.normal(self.absorb_tag, self.dimensions)
 
-        from grudge.symbolic import RestrictToBoundary, Field
 
-        e, h = self.split_eh(self.field_placeholder(w))
+        e, h = self.split_eh(w)
 
         if self.fixed_material:
             epsilon = self.epsilon
             mu = self.mu
-        else:
-            epsilon = cse(
-                    RestrictToBoundary(self.absorb_tag)(Field("epsilon")))
-            mu = cse(
-                    RestrictToBoundary(self.absorb_tag)(Field("mu")))
 
         absorb_Z = (mu/epsilon)**0.5
         absorb_Y = 1/absorb_Z
 
-        absorb_e = RestrictToBoundary(self.absorb_tag)(e)
-        absorb_h = RestrictToBoundary(self.absorb_tag)(h)
+        absorb_e = sym.cse(sym.interp("vol", self.absorb_tag)(e))
+        absorb_h = sym.cse(sym.interp("vol", self.absorb_tag)(h))
 
         bc = join_fields(
                 absorb_e + 1/2*(self.space_cross_h(absorb_normal, self.space_cross_e(
@@ -281,17 +238,12 @@ class MaxwellOperator(HyperbolicOperator):
 
         return bc
 
-    def incident_bc(self, w=None):
+    def incident_bc(self, w):
         "Flux terms for incident boundary conditions"
         # NOTE: Untested for inhomogeneous materials, but would usually be
         # physically meaningless anyway (are there exceptions to this?)
 
-        e, h = self.split_eh(self.field_placeholder(w))
-        if not self.fixed_material:
-            from warnings import warn
-            if self.incident_tag != BTAG_NONE:
-                warn("Incident boundary conditions assume homogeneous"
-                     " background material, results may be unphysical")
+        e, h = self.split_eh(w)
 
         from grudge.tools import count_subset
         fld_cnt = count_subset(self.get_eh_subset())
@@ -301,7 +253,7 @@ class MaxwellOperator(HyperbolicOperator):
         if is_zero(incident_bc_data):
             return make_obj_array([0]*fld_cnt)
         else:
-            return cse(-incident_bc_data)
+            return sym.cse(-incident_bc_data)
 
     def sym_operator(self, w=None):
         """The full operator template - the high level description of
@@ -310,23 +262,9 @@ class MaxwellOperator(HyperbolicOperator):
         Combines the relevant operator templates for spatial
         derivatives, flux, boundary conditions etc.
         """
-        w = self.field_placeholder(w)
-
-        if self.fixed_material:
-            flux_w = w
-        else:
-            epsilon = self.epsilon
-            mu = self.mu
-
-            flux_w = join_fields(epsilon, mu, w)
-
-        from grudge.symbolic import BoundaryPair, \
-                InverseMassOperator, get_flux_operator
-
-        flux_op = get_flux_operator(self.flux(self.flux_type))
-        bdry_flux_op = get_flux_operator(self.flux(self.bdry_flux_type))
+        from grudge.tools import count_subset
+        w = sym.make_sym_array("w", count_subset(self.get_eh_subset()))
 
-        from grudge.tools.indexing import count_subset
         elec_components = count_subset(self.get_eh_subset()[0:3])
         mag_components = count_subset(self.get_eh_subset()[3:6])
 
@@ -334,10 +272,6 @@ class MaxwellOperator(HyperbolicOperator):
             # need to check this
             material_divisor = (
                     [self.epsilon]*elec_components+[self.mu]*mag_components)
-        else:
-            material_divisor = join_fields(
-                    [epsilon]*elec_components,
-                    [mu]*mag_components)
 
         tags_and_bcs = [
                 (self.pec_tag, self.pec_bc(w)),
@@ -346,25 +280,19 @@ class MaxwellOperator(HyperbolicOperator):
                 (self.incident_tag, self.incident_bc(w)),
                 ]
 
-        def make_flux_bc_vector(tag, bc):
-            if self.fixed_material:
-                return bc
-            else:
-                from grudge.symbolic import RestrictToBoundary
-                return join_fields(
-                        cse(RestrictToBoundary(tag)(epsilon)),
-                        cse(RestrictToBoundary(tag)(mu)),
-                        bc)
+
+        def flux(pair):
+            return sym.interp(pair.dd, "all_faces")(self.flux(pair))
 
         return (
                 - self.local_derivatives(w)
-                + InverseMassOperator()(
-                    flux_op(flux_w)
-                    + sum(
-                        bdry_flux_op(BoundaryPair(
-                            flux_w, make_flux_bc_vector(tag, bc), tag))
-                        for tag, bc in tags_and_bcs))
-                    ) / material_divisor
+                + sym.InverseMassOperator()(sym.FaceMassOperator()(
+                    flux(sym.int_tpair(w))
+                    #+ sum(
+                        #flux(sym.bv_tpair(tag, w, bc))
+                            #for tag, bc in tags_and_bcs)
+                    ))) / material_divisor
+
 
     def bind(self, discr):
         "Convert the abstract operator template into compiled code."
@@ -445,11 +373,8 @@ class MaxwellOperator(HyperbolicOperator):
         e_idx, h_idx = self.partial_to_eh_subsets()
         e, h = w[e_idx], w[h_idx]
 
-        from grudge.flux import FluxVectorPlaceholder as FVP
-        if isinstance(w, FVP):
-            return FVP(scalars=e), FVP(scalars=h)
-        else:
-            return make_obj_array(e), make_obj_array(h)
+        return e,h
+        #return make_obj_array(e), make_obj_array(h)
 
     def get_eh_subset(self):
         """Return a 6-tuple of :class:`bool` objects indicating whether field
@@ -476,6 +401,14 @@ class MaxwellOperator(HyperbolicOperator):
             raise ValueError("max_eigenvalue is no longer supported for "
                     "variable-coefficient problems--use max_eigenvalue_expr")
 
+    def check_bc_coverage(self, mesh):
+        from meshmode.mesh import check_bc_coverage
+        check_bc_coverage(mesh, [
+            self.pec_tag,
+            self.pmc_tag,
+            self.absorb_tag,
+            self.incident_tag])
+
 
 class TMMaxwellOperator(MaxwellOperator):
     """A 2D TM Maxwell operator with PEC boundaries.

grudge/models/wave.py

@@ -180,187 +180,303 @@ 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))))
+
+    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)
 
-        return flux
+        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")
 
-    def bind_characteristic_velocity(self, discr):
-        from grudge.symbolic.operators import ElementwiseMaxOperator
+        # radiation BCs -------------------------------------------------------
+        rad_normal = sym.normal(self.radiation_tag, d)
 
-        compiled = discr.compile(ElementwiseMaxOperator()(self.c))
+        rad_u = sym.cse(sym.interp("vol", self.radiation_tag)(u))
+        rad_v = sym.cse(sym.interp("vol", self.radiation_tag)(v))
 
-        def do(t, w, **extra_context):
-            return compiled(t=t, w=w, **extra_context)
+        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")
 
-        return do
+        # entire operator -----------------------------------------------------
+        nabla = sym.nabla(d)
 
-    def sym_operator(self, with_sensor=False):
-        d = self.ambient_dim
+        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)
+                    )
 
-        w = sym.make_sym_array("w", d+1)
-        u = w[0]
-        v = w[1:]
 
-        flux_w = join_fields(self.c, 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
+
+        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)
+
 
-        # {{{ boundary conditions
+# }}}
+
+
+# {{{ variable-velocity
+
+class VariableCoefficientWeakWaveOperator(HyperbolicOperator):
+    """This operator discretizes the wave equation
+    :math:`\\partial_t^2 u = c^2 \\Delta u`.
 
-        # 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
+    To be precise, we discretize the hyperbolic system
 
-        dir_bc = join_fields(dir_c, -dir_u, dir_v)
+    .. math::
 
-        # 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
+        \partial_t u - c \\nabla \\cdot v = 0
 
-        neu_bc = join_fields(neu_c, neu_u, -neu_v)
+        \partial_t v - c \\nabla u = 0
 
-        # Radiation
-        rad_normal = sym.make_normal(self.radiation_tag, d)
+    The sign of :math:`v` determines whether we discretize the forward or the
+    backward wave equation.
 
-        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
+    :math:`c` is assumed to be constant across all space.
+    """
 
-        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)
-                )
+    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
+
+        self.sign = sym.If(sym.Comparison(
+                            self.c, ">", 0),
+                                            1, -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)
+
+        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
+        if 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))))
 
-        # {{{ diffusion -------------------------------------------------------
-        from pytools.obj_array import with_object_array_or_scalar
+    def sym_operator(self):
+        d = self.ambient_dim
 
-        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
+        w = sym.make_sym_array("w", d+1)
+        u = w[0]
+        v = w[1:]
 
-                if with_sensor:
-                    diffusion_coeff += sym.Field("sensor")
+        # boundary conditions -------------------------------------------------
 
-                from grudge.second_order import SecondDerivativeTarget
+        # 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.")
 
-                # strong_form here allows the reuse the value of grad u.
-                grad_tgt = SecondDerivativeTarget(
-                        self.ambient_dim, strong_form=True,
-                        operand=arg)
+            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)
 
-                self.diffusion_scheme.grad(grad_tgt, bc_getter=None,
-                        dirichlet_tags=[], neumann_tags=[])
+        dir_bc = sym.cse(dir_bc, "dir_bc")
 
-                div_tgt = SecondDerivativeTarget(
-                        self.ambient_dim, strong_form=False,
-                        operand=diffusion_coeff*grad_tgt.minv_all)
+        # 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")
 
-                self.diffusion_scheme.div(div_tgt,
-                        bc_getter=None,
-                        dirichlet_tags=[], neumann_tags=[])
+        # radiation BCs -------------------------------------------------------
+        rad_normal = sym.normal(self.radiation_tag, d)
 
-                return div_tgt.minv_all
-            else:
-                return 0
+        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 +485,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)
+
 
 # }}}
 

grudge/tools.py

@@ -36,3 +36,105 @@ def is_zero(x):
         return x == 0
     else:
         return False
+
+def levi_civita(tuple):
+    """Compute an entry of the Levi-Civita tensor for the indices *tuple*.
+    Only three-tuples are supported for now.
+    """
+    if len(tuple) == 3:
+        if tuple in [(0,1,2), (2,0,1), (1,2,0)]:
+            return 1
+        elif tuple in [(2,1,0), (0,2,1), (1,0,2)]:
+            return -1
+        else:
+            return 0
+    else:
+        raise NotImplementedError
+
+
+class SubsettableCrossProduct:
+    """A cross product that can operate on an arbitrary subsets of its
+    two operands and return an arbitrary subset of its result.
+    """
+
+    full_subset = (True, True, True)
+
+    def __init__(self, op1_subset=full_subset, op2_subset=full_subset, result_subset=full_subset):
+        """Construct a subset-able cross product.
+        :param op1_subset: The subset of indices of operand 1 to be taken into account.
+          Given as a 3-sequence of bools.
+        :param op2_subset: The subset of indices of operand 2 to be taken into account.
+          Given as a 3-sequence of bools.
+        :param result_subset: The subset of indices of the result that are calculated.
+          Given as a 3-sequence of bools.
+        """
+        def subset_indices(subset):
+            return [i for i, use_component in enumerate(subset)
+                    if use_component]
+
+        self.op1_subset = op1_subset
+        self.op2_subset = op2_subset
+        self.result_subset = result_subset
+
+        import pymbolic
+        op1 = pymbolic.var("x")
+        op2 = pymbolic.var("y")
+
+        self.functions = []
+        self.component_lcjk = []
+        for i, use_component in enumerate(result_subset):
+            if use_component:
+                this_expr = 0
+                this_component = []
+                for j, j_real in enumerate(subset_indices(op1_subset)):
+                    for k, k_real in enumerate(subset_indices(op2_subset)):
+                        lc = levi_civita((i, j_real, k_real))
+                        if lc != 0:
+                            this_expr += lc*op1.index(j)*op2.index(k)
+                            this_component.append((lc, j, k))
+                self.functions.append(pymbolic.compile(this_expr,
+                    variables=[op1, op2]))
+                self.component_lcjk.append(this_component)
+
+    def __call__(self, x, y, three_mult=None):
+        """Compute the subsetted cross product on the indexables *x* and *y*.
+        :param three_mult: a function of three arguments *sign, xj, yk*
+          used in place of the product *sign*xj*yk*. Defaults to just this
+          product if not given.
+        """
+        from pytools.obj_array import join_fields
+        if three_mult is None:
+            return join_fields(*[f(x, y) for f in self.functions])
+        else:
+            return join_fields(
+                    *[sum(three_mult(lc, x[j], y[k]) for lc, j, k in lcjk)
+                    for lcjk in self.component_lcjk])
+
+
+
+
+cross = SubsettableCrossProduct()
+
+def count_subset(subset):
+    from pytools import len_iterable
+    return len_iterable(uc for uc in subset if uc)
+
+def partial_to_all_subset_indices(subsets, base=0):
+    """Takes a sequence of bools and generates it into an array of indices
+    to be used to insert the subset into the full set.
+    Example:
+    >>> list(partial_to_all_subset_indices([[False, True, True], [True,False,True]]))
+    [array([0 1]), array([2 3]
+    """
+
+    idx = base
+    for subset in subsets:
+        result = []
+        for is_in in subset:
+            if is_in:
+                result.append(idx)
+                idx += 1
+
+        yield np.array(result, dtype=np.intp)
+
+