diff --git a/examples/wave/var-propagation-speed.py b/examples/wave/var-propagation-speed.py
index 3f4b9927c0436896d94e95e78f2d1325ec4a2045..24719059263217f760745f84afa2851668775450 100644
--- a/examples/wave/var-propagation-speed.py
+++ b/examples/wave/var-propagation-speed.py
@@ -1,11 +1,8 @@
-"""Variable-coefficient wave propagation."""
+"""Minimal example of a grudge driver."""
-from __future__ import division
-from __future__ import absolute_import
-from __future__ import print_function
-from six.moves import range
+from __future__ import division, print_function
-__copyright__ = "Copyright (C) 2009 Andreas Kloeckner"
+__copyright__ = "Copyright (C) 2015 Andreas Kloeckner"
__license__ = """
Permission is hereby granted, free of charge, to any person obtaining a copy
@@ -29,171 +26,105 @@ THE SOFTWARE.
import numpy as np
-from grudge.mesh import BTAG_ALL, BTAG_NONE
-
-
-def main(write_output=True,
- dir_tag=BTAG_NONE,
- neu_tag=BTAG_NONE,
- rad_tag=BTAG_ALL,
- flux_type_arg="upwind"):
- from math import sin, cos, pi, exp, sqrt # noqa
-
- from grudge.backends import guess_run_context
- rcon = guess_run_context()
-
- dim = 2
-
- if dim == 1:
- if rcon.is_head_rank:
- from grudge.mesh.generator import make_uniform_1d_mesh
- mesh = make_uniform_1d_mesh(-10, 10, 500)
- elif dim == 2:
- from grudge.mesh.generator import make_rect_mesh
- if rcon.is_head_rank:
- mesh = make_rect_mesh(a=(-1, -1), b=(1, 1), max_area=0.003)
- elif dim == 3:
- if rcon.is_head_rank:
- from grudge.mesh.generator import make_ball_mesh
- mesh = make_ball_mesh(max_volume=0.0005)
- else:
- raise RuntimeError("bad number of dimensions")
-
- if rcon.is_head_rank:
- print("%d elements" % len(mesh.elements))
- mesh_data = rcon.distribute_mesh(mesh)
- else:
- mesh_data = rcon.receive_mesh()
-
- discr = rcon.make_discretization(mesh_data, order=4)
-
- from grudge.timestep.runge_kutta import LSRK4TimeStepper
- stepper = LSRK4TimeStepper()
-
- from grudge.visualization import VtkVisualizer
- if write_output:
- vis = VtkVisualizer(discr, rcon, "fld")
-
- source_center = np.array([0.7, 0.4])
- source_width = 1/16
- source_omega = 3
-
- import grudge.symbolic as sym
- sym_x = sym.nodes(2)
- sym_source_center_dist = sym_x - source_center
+import pyopencl as cl
+from grudge.shortcuts import set_up_rk4
+from grudge import sym, bind, Discretization
- from grudge.models.wave import VariableVelocityStrongWaveOperator
- op = VariableVelocityStrongWaveOperator(
- c=sym.If(sym.Comparison(
- np.dot(sym_x, sym_x), "<", 0.4**2),
- 1, 0.5),
- dimensions=discr.dimensions,
- source=
- sym.CFunction("sin")(source_omega*sym.ScalarParameter("t"))
- * sym.CFunction("exp")(
- -np.dot(sym_source_center_dist, sym_source_center_dist)
- / source_width**2),
- dirichlet_tag=dir_tag,
- neumann_tag=neu_tag,
- radiation_tag=rad_tag,
- flux_type=flux_type_arg
- )
-
- from grudge.tools import join_fields
- fields = join_fields(discr.volume_zeros(),
- [discr.volume_zeros() for i in range(discr.dimensions)])
-
- # {{{ diagnostics setup
-
- from pytools.log import LogManager, \
- add_general_quantities, \
- add_simulation_quantities, \
- add_run_info
-
- if write_output:
- log_file_name = "wave.dat"
- 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)
-
- from pytools.log import IntervalTimer
- vis_timer = IntervalTimer("t_vis", "Time spent visualizing")
- logmgr.add_quantity(vis_timer)
- stepper.add_instrumentation(logmgr)
-
- from grudge.log import LpNorm
- u_getter = lambda: fields[0]
- logmgr.add_quantity(LpNorm(u_getter, discr, 1, name="l1_u"))
- logmgr.add_quantity(LpNorm(u_getter, discr, name="l2_u"))
-
- logmgr.add_watches(["step.max", "t_sim.max", "l2_u", "t_step.max"])
-
- # }}}
-
- # {{{ timestep loop
-
- rhs = op.bind(discr)
- try:
- from grudge.timestep.stability import \
- approximate_rk4_relative_imag_stability_region
- max_dt = (
- 1/discr.compile(op.max_eigenvalue_expr())()
- * discr.dt_non_geometric_factor()
- * discr.dt_geometric_factor()
- * approximate_rk4_relative_imag_stability_region(stepper))
- if flux_type_arg == "central":
- max_dt *= 0.25
-
- from grudge.timestep import times_and_steps
- step_it = times_and_steps(final_time=3, logmgr=logmgr,
- max_dt_getter=lambda t: max_dt)
-
- for step, t, dt in step_it:
- if step % 10 == 0 and write_output:
- visf = vis.make_file("fld-%04d" % step)
-
- vis.add_data(visf,
- [
- ("u", fields[0]),
- ("v", fields[1:]),
- ],
- time=t,
- step=step)
- visf.close()
-
- fields = stepper(fields, t, dt, rhs)
- assert discr.norm(fields) < 1
- finally:
- if write_output:
- vis.close()
+def main(write_output=True, order=4):
+ cl_ctx = cl.create_some_context()
+ queue = cl.CommandQueue(cl_ctx)
- logmgr.close()
- discr.close()
+ 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)
- # }}}
-
-if __name__ == "__main__":
- main(flux_type_arg="upwind")
+ 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
-# entry points for py.test ----------------------------------------------------
-def test_var_velocity_wave():
- from pytools.test import mark_test
- mark_long = mark_test.long
+ 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()
- for flux_type in ["upwind", "central"]:
- yield ("dirichlet var-v wave equation with %s flux" % flux_type,
- mark_long(main),
- False, BTAG_ALL, BTAG_NONE, TAG_NONE, flux_type)
- yield ("neumann var-v wave equation", mark_long(main),
- False, BTAG_NONE, BTAG_ALL, TAG_NONE)
- yield ("radiation-bc var-v wave equation", mark_long(main),
- False, BTAG_NONE, TAG_NONE, BTAG_ALL)
-# vim: foldmethod=marker
+if __name__ == "__main__":
+ main()
diff --git a/grudge/execution.py b/grudge/execution.py
index 0241c82b532c9167c60362917a08c73123c385df..f78c345f3c87f4c495afa3685b15680c06fae0cd 100644
--- a/grudge/execution.py
+++ b/grudge/execution.py
@@ -119,12 +119,37 @@ class ExecutionMapper(mappers.Evaluator,
then = self.rec(expr.then)
else_ = self.rec(expr.else_)
- result = cl.array.empty_like(then, queue=self.queue,
- allocator=self.bound_op.allocator)
- cl.array.if_positive(bool_crit, then, else_, out=result,
- queue=self.queue)
+ import pymbolic.primitives as p
- return result
+ var = p.Variable
+
+ i = var("i") #sym.Variable("i")
+ if isinstance(then, pyopencl.array.Array):
+ sym_then = var("a")[i] #sym.Variable("a")[i]
+ else:
+ sym_then = var("a") # sym.Variable("a")
+ then = np.float64(then)
+
+ if isinstance(else_, pyopencl.array.Array):
+ sym_else = var("b")[i] # sym.Variable("b")[i]
+ else:
+ sym_else = var("b") # sym.Variable("b")
+ else_ = np.float64(else_)
+
+ @memoize_in(self.bound_op, "map_if_knl")
+ def knl():
+ knl = lp.make_kernel(
+ "{[i]: 0<=i<n}",
+ [
+ lp.Assignment(var("out")[i], p.If(var("crit")[i], sym_then, sym_else))
+
+ #lp.Assignment(sym.Variable("out")[i], sym.If(sym.Variable("crit")[i], sym_then, sym_else))
+ ])
+ return lp.split_iname(knl, "i", 128, outer_tag="g.0", inner_tag="l.0")
+
+ evt, (out,) = knl()(self.queue, crit=bool_crit, a=then, b=else_)
+
+ return out
def map_ref_diff_base(self, op, field_expr):
raise NotImplementedError(
diff --git a/grudge/models/wave.py b/grudge/models/wave.py
index 98f0b4407324b82895da2040503de692849f9743..6ea1d923287761bba5efbb4998beda4a519e6318 100644
--- a/grudge/models/wave.py
+++ b/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)
+
# }}}