diff --git a/examples/wave/wave-eager-mpi.py b/examples/wave/wave-eager-mpi.py
index 59de0ff90937c634da0a4d3c28b2f3dc031b675b..4b408bb599c5002a32325f5befd9ab66ab02cee9 100644
--- a/examples/wave/wave-eager-mpi.py
+++ b/examples/wave/wave-eager-mpi.py
@@ -43,26 +43,26 @@ from mpi4py import MPI
# {{{ wave equation bits
+def scalar(arg):
+ return make_obj_array([arg])
+
+
def wave_flux(discr, c, w_tpair):
u = w_tpair[0]
v = w_tpair[1:]
normal = thaw(u.int.array_context, discr.normal(w_tpair.dd))
- def normal_times(scalar):
- # workaround for object array behavior
- return make_obj_array([ni*scalar for ni in normal])
-
flux_weak = flat_obj_array(
np.dot(v.avg, normal),
- normal_times(u.avg),
+ normal*scalar(u.avg),
)
# upwind
v_jump = np.dot(normal, v.int-v.ext)
flux_weak -= flat_obj_array(
0.5*(u.int-u.ext),
- 0.5*normal_times(v_jump),
+ 0.5*normal*scalar(v_jump),
)
return discr.project(w_tpair.dd, "all_faces", c*flux_weak)
@@ -175,7 +175,7 @@ def main():
[discr.zeros(actx) for i in range(discr.dim)]
)
- vis = make_visualizer(discr, discr.order+3 if dim == 2 else discr.order)
+ vis = make_visualizer(discr, order+3 if dim == 2 else order)
def rhs(t, w):
return wave_operator(discr, c=1, w=w)
diff --git a/examples/wave/wave-eager-var-velocity.py b/examples/wave/wave-eager-var-velocity.py
new file mode 100644
index 0000000000000000000000000000000000000000..721211314dfd5931374f0347e90435c624baf129
--- /dev/null
+++ b/examples/wave/wave-eager-var-velocity.py
@@ -0,0 +1,211 @@
+from __future__ import division, print_function
+
+__copyright__ = "Copyright (C) 2020 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 numpy.linalg as la # noqa
+import pyopencl as cl
+
+from pytools.obj_array import flat_obj_array, make_obj_array
+
+from meshmode.array_context import PyOpenCLArrayContext
+from meshmode.dof_array import thaw
+
+from meshmode.mesh import BTAG_ALL, BTAG_NONE # noqa
+
+from grudge.eager import EagerDGDiscretization, interior_trace_pair
+from grudge.shortcuts import make_visualizer
+from grudge.symbolic.primitives import TracePair, QTAG_NONE, DOFDesc
+
+
+# {{{ wave equation bits
+
+def scalar(arg):
+ return make_obj_array([arg])
+
+
+def wave_flux(discr, c, w_tpair):
+ dd = w_tpair.dd
+ dd_quad = dd.with_qtag("vel_prod")
+
+ u = w_tpair[0]
+ v = w_tpair[1:]
+
+ normal = thaw(u.int.array_context, discr.normal(dd))
+
+ flux_weak = flat_obj_array(
+ np.dot(v.avg, normal),
+ normal*scalar(u.avg),
+ )
+
+ # upwind
+ v_jump = np.dot(normal, v.int-v.ext)
+ flux_weak -= flat_obj_array(
+ 0.5*(u.int-u.ext),
+ 0.5*normal*scalar(v_jump),
+ )
+
+ # FIMXE this flux is only correct for continuous c
+ dd_allfaces_quad = dd_quad.with_dtag("all_faces")
+ c_quad = discr.project("vol", dd_quad, c)
+ flux_quad = discr.project(dd, dd_quad, flux_weak)
+
+ return discr.project(dd_quad, dd_allfaces_quad, scalar(c_quad)*flux_quad)
+
+
+def wave_operator(discr, c, w):
+ u = w[0]
+ v = w[1:]
+
+ dir_u = discr.project("vol", BTAG_ALL, u)
+ dir_v = discr.project("vol", BTAG_ALL, v)
+ dir_bval = flat_obj_array(dir_u, dir_v)
+ dir_bc = flat_obj_array(-dir_u, dir_v)
+
+ dd_quad = DOFDesc("vol", "vel_prod")
+ c_quad = discr.project("vol", dd_quad, c)
+ w_quad = discr.project("vol", dd_quad, w)
+ u_quad = w_quad[0]
+ v_quad = w_quad[1:]
+
+ dd_allfaces_quad = DOFDesc("all_faces", "vel_prod")
+
+ # FIXME Fix sign issue
+ return (
+ discr.inverse_mass(
+ flat_obj_array(
+ discr.weak_div(dd_quad, scalar(c_quad)*v_quad),
+ discr.weak_grad(dd_quad, c_quad*u_quad)
+ )
+ - # noqa: W504
+ discr.face_mass(
+ dd_allfaces_quad,
+ wave_flux(discr, c=c, w_tpair=interior_trace_pair(discr, w))
+ + wave_flux(discr, c=c, w_tpair=TracePair(
+ BTAG_ALL, dir_bval, dir_bc))
+ ))
+ )
+
+# }}}
+
+
+def rk4_step(y, t, h, f):
+ k1 = f(t, y)
+ k2 = f(t+h/2, y + h/2*k1)
+ k3 = f(t+h/2, y + h/2*k2)
+ k4 = f(t+h, y + h*k3)
+ return y + h/6*(k1 + 2*k2 + 2*k3 + k4)
+
+
+def bump(actx, discr, t=0, width=0.05, center=None):
+ if center is None:
+ center = np.array([0.2, 0.35, 0.1])
+
+ center = center[:discr.dim]
+ source_omega = 3
+
+ nodes = thaw(actx, discr.nodes())
+ center_dist = flat_obj_array([
+ nodes[i] - center[i]
+ for i in range(discr.dim)
+ ])
+
+ return (
+ np.cos(source_omega*t)
+ * actx.np.exp(
+ -np.dot(center_dist, center_dist)
+ / width**2))
+
+
+def main():
+ cl_ctx = cl.create_some_context()
+ queue = cl.CommandQueue(cl_ctx)
+ actx = PyOpenCLArrayContext(queue)
+
+ dim = 2
+ nel_1d = 16
+ from meshmode.mesh.generation import generate_regular_rect_mesh
+ mesh = generate_regular_rect_mesh(
+ a=(-0.5,)*dim,
+ b=(0.5,)*dim,
+ n=(nel_1d,)*dim)
+
+ order = 3
+
+ if dim == 2:
+ # no deep meaning here, just a fudge factor
+ dt = 0.75/(nel_1d*order**2)
+ elif dim == 3:
+ # no deep meaning here, just a fudge factor
+ dt = 0.45/(nel_1d*order**2)
+ else:
+ raise ValueError("don't have a stable time step guesstimate")
+
+ print("%d elements" % mesh.nelements)
+
+ from meshmode.discretization.poly_element import \
+ QuadratureSimplexGroupFactory, \
+ PolynomialWarpAndBlendGroupFactory
+ discr = EagerDGDiscretization(actx, mesh,
+ quad_tag_to_group_factory={
+ QTAG_NONE: PolynomialWarpAndBlendGroupFactory(order),
+ "vel_prod": QuadratureSimplexGroupFactory(3*order),
+ })
+
+ # bounded above by 1
+ c = 0.2 + 0.8*bump(actx, discr, center=np.zeros(3), width=0.5)
+
+ fields = flat_obj_array(
+ bump(actx, discr, ),
+ [discr.zeros(actx) for i in range(discr.dim)]
+ )
+
+ vis = make_visualizer(discr, order+3 if dim == 2 else order)
+
+ def rhs(t, w):
+ return wave_operator(discr, c=c, w=w)
+
+ t = 0
+ t_final = 3
+ istep = 0
+ while t < t_final:
+ fields = rk4_step(fields, t, dt, rhs)
+
+ if istep % 10 == 0:
+ print(istep, t, discr.norm(fields[0]))
+ vis.write_vtk_file("fld-wave-eager-var-velocity-%04d.vtu" % istep,
+ [
+ ("c", c),
+ ("u", fields[0]),
+ ("v", fields[1:]),
+ ])
+
+ t += dt
+ istep += 1
+
+
+if __name__ == "__main__":
+ main()
+
+# vim: foldmethod=marker
diff --git a/examples/wave/wave-eager.py b/examples/wave/wave-eager.py
index ba1e4af7ff77dc214c7c03207415e16c09164ab6..e78214313bfb8ea0ab9e547bbc37d6ebc0aaaa5f 100644
--- a/examples/wave/wave-eager.py
+++ b/examples/wave/wave-eager.py
@@ -41,26 +41,26 @@ from grudge.symbolic.primitives import TracePair
# {{{ wave equation bits
+def scalar(arg):
+ return make_obj_array([arg])
+
+
def wave_flux(discr, c, w_tpair):
u = w_tpair[0]
v = w_tpair[1:]
normal = thaw(u.int.array_context, discr.normal(w_tpair.dd))
- def normal_times(scalar):
- # workaround for object array behavior
- return make_obj_array([ni*scalar for ni in normal])
-
flux_weak = flat_obj_array(
np.dot(v.avg, normal),
- normal_times(u.avg),
+ normal*scalar(u.avg),
)
# upwind
v_jump = np.dot(normal, v.int-v.ext)
flux_weak -= flat_obj_array(
0.5*(u.int-u.ext),
- 0.5*normal_times(v_jump),
+ 0.5*normal*scalar(v_jump),
)
return discr.project(w_tpair.dd, "all_faces", c*flux_weak)
@@ -151,7 +151,7 @@ def main():
[discr.zeros(actx) for i in range(discr.dim)]
)
- vis = make_visualizer(discr, discr.order+3 if dim == 2 else discr.order)
+ vis = make_visualizer(discr, order+3 if dim == 2 else order)
def rhs(t, w):
return wave_operator(discr, c=1, w=w)
diff --git a/grudge/discretization.py b/grudge/discretization.py
index a8bee7e11a046aaeaf810f3f47c0e01ae6df995d..7c0573c2e5ec3515843b6f8c377ad10301256aaa 100644
--- a/grudge/discretization.py
+++ b/grudge/discretization.py
@@ -48,8 +48,8 @@ class DGDiscretizationWithBoundaries(DiscretizationBase):
.. automethod :: zeros
"""
- def __init__(self, array_context, mesh, order, quad_tag_to_group_factory=None,
- mpi_communicator=None):
+ def __init__(self, array_context, mesh, order=None,
+ quad_tag_to_group_factory=None, mpi_communicator=None):
"""
:param quad_tag_to_group_factory: A mapping from quadrature tags (typically
strings--but may be any hashable/comparable object) to a
@@ -61,10 +61,27 @@ class DGDiscretizationWithBoundaries(DiscretizationBase):
self._setup_actx = array_context
+ from meshmode.discretization.poly_element import \
+ PolynomialWarpAndBlendGroupFactory
+
if quad_tag_to_group_factory is None:
- quad_tag_to_group_factory = {}
+ if order is None:
+ raise TypeError("one of 'order' and "
+ "'quad_tag_to_group_factory' must be given")
+
+ quad_tag_to_group_factory = {
+ sym.QTAG_NONE: PolynomialWarpAndBlendGroupFactory(order=order)}
+ else:
+ if order is not None:
+ quad_tag_to_group_factory = quad_tag_to_group_factory.copy()
+ if sym.QTAG_NONE in quad_tag_to_group_factory:
+ raise ValueError("if 'order' is given, "
+ "'quad_tag_to_group_factory' must not have a "
+ "key of QTAG_NONE")
+
+ quad_tag_to_group_factory[sym.QTAG_NONE] = \
+ PolynomialWarpAndBlendGroupFactory(order=order)
- self.order = order
self.quad_tag_to_group_factory = quad_tag_to_group_factory
from meshmode.discretization import Discretization
@@ -268,13 +285,7 @@ class DGDiscretizationWithBoundaries(DiscretizationBase):
if quadrature_tag is None:
quadrature_tag = sym.QTAG_NONE
- from meshmode.discretization.poly_element import \
- PolynomialWarpAndBlendGroupFactory
-
- if quadrature_tag is not sym.QTAG_NONE:
- return self.quad_tag_to_group_factory[quadrature_tag]
- else:
- return PolynomialWarpAndBlendGroupFactory(order=self.order)
+ return self.quad_tag_to_group_factory[quadrature_tag]
@memoize_method
def _quad_volume_discr(self, quadrature_tag):
@@ -375,5 +386,16 @@ class DGDiscretizationWithBoundaries(DiscretizationBase):
where is None
or where == sym.VTAG_ALL)
+ @property
+ def order(self):
+ from warnings import warn
+ warn("DGDiscretizationWithBoundaries.order is deprecated, "
+ "consider the orders of element groups instead. "
+ "'order' will go away in 2021.",
+ DeprecationWarning, stacklevel=2)
+
+ from pytools import single_valued
+ return single_valued(egrp.order for egrp in self._volume_discr.groups)
+
# vim: foldmethod=marker
diff --git a/grudge/eager.py b/grudge/eager.py
index d2525a7d3418872bf189f6021cad7ca8fd25b63b..ffea3b599312f30d82b300698e2dd207d86cbec7 100644
--- a/grudge/eager.py
+++ b/grudge/eager.py
@@ -68,16 +68,33 @@ class EagerDGDiscretization(DGDiscretizationWithBoundaries):
self.grad(vec_i)[i] for i, vec_i in enumerate(vecs))
@memoize_method
- def _bound_weak_grad(self):
- return bind(self, sym.stiffness_t(self.dim) * sym.Variable("u"),
+ def _bound_weak_grad(self, dd):
+ return bind(self,
+ sym.stiffness_t(self.dim, dd_in=dd) * sym.Variable("u", dd),
local_only=True)
- def weak_grad(self, vec):
- return self._bound_weak_grad()(u=vec)
+ def weak_grad(self, *args):
+ if len(args) == 1:
+ vec, = args
+ dd = sym.DOFDesc("vol", sym.QTAG_NONE)
+ elif len(args) == 2:
+ dd, vec = args
+ else:
+ raise TypeError("invalid number of arguments")
+
+ return self._bound_weak_grad(dd)(u=vec)
+
+ def weak_div(self, *args):
+ if len(args) == 1:
+ vecs, = args
+ dd = sym.DOFDesc("vol", sym.QTAG_NONE)
+ elif len(args) == 2:
+ dd, vecs = args
+ else:
+ raise TypeError("invalid number of arguments")
- def weak_div(self, vecs):
return sum(
- self.weak_grad(vec_i)[i] for i, vec_i in enumerate(vecs))
+ self.weak_grad(dd, vec_i)[i] for i, vec_i in enumerate(vecs))
@memoize_method
def normal(self, dd):
@@ -104,18 +121,26 @@ class EagerDGDiscretization(DGDiscretizationWithBoundaries):
return self._bound_inverse_mass()(u=vec)
@memoize_method
- def _bound_face_mass(self):
- u = sym.Variable("u", dd=sym.as_dofdesc("all_faces"))
- return bind(self, sym.FaceMassOperator()(u), local_only=True)
+ def _bound_face_mass(self, dd):
+ u = sym.Variable("u", dd=dd)
+ return bind(self, sym.FaceMassOperator(dd_in=dd)(u), local_only=True)
+
+ def face_mass(self, *args):
+ if len(args) == 1:
+ vec, = args
+ dd = sym.DOFDesc("all_faces", sym.QTAG_NONE)
+ elif len(args) == 2:
+ dd, vec = args
+ else:
+ raise TypeError("invalid number of arguments")
- def face_mass(self, vec):
if (isinstance(vec, np.ndarray)
and vec.dtype.char == "O"
and not isinstance(vec, DOFArray)):
return obj_array_vectorize(
- lambda el: self.face_mass(el), vec)
+ lambda el: self.face_mass(dd, el), vec)
- return self._bound_face_mass()(u=vec)
+ return self._bound_face_mass(dd)(u=vec)
@memoize_method
def _norm(self, p):