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):