diff --git a/doc/conf.py b/doc/conf.py index 866653625acb09a79319be6bec7c07872081f38f..5d9f499e4c3a0cc10b43ec14611fe88365bde549 100644 --- a/doc/conf.py +++ b/doc/conf.py @@ -280,5 +280,6 @@ intersphinx_mapping = { 'https://documen.tician.de/pyopencl': None, 'https://documen.tician.de/meshpy': None, 'https://documen.tician.de/modepy': None, - 'https://documen.tician.de/loopy': None + 'https://documen.tician.de/loopy': None, + 'https://tisaac.gitlab.io/recursivenodes/': None, } diff --git a/meshmode/array_context.py b/meshmode/array_context.py index 6118e9d5530c1ee46aff3ff75ed4d0413698420b..88d3ecf412f033833272d696de2e70c9b35266f6 100644 --- a/meshmode/array_context.py +++ b/meshmode/array_context.py @@ -153,6 +153,14 @@ class ArrayContext: @memoize_method def _get_scalar_func_loopy_program(self, name, nargs, naxes): + if name == "arctan2": + name = "atan2" + elif name == "atan2": + from warnings import warn + warn("'atan2' in ArrayContext.np is deprecated. Use 'arctan2', " + "as in numpy2. This will be disallowed in 2021.", + DeprecationWarning, stacklevel=3) + from pymbolic import var var_names = ["i%d" % i for i in range(naxes)] @@ -216,6 +224,11 @@ class _PyOpenCLFakeNumpyNamespace(_BaseFakeNumpyNamespace): return super().__getattr__(name) + @obj_array_vectorized_n_args + def where(self, criterion, then, else_): + import pyopencl.array as cl_array + return cl_array.if_positive(criterion != 0, then, else_) + class PyOpenCLArrayContext(ArrayContext): """ diff --git a/meshmode/discretization/poly_element.py b/meshmode/discretization/poly_element.py index 5de2102cb4c227329555e2d7cbb1d24ad3f5b7b0..55cce2eab1c9b994225b5ea8caa91dfa21912523 100644 --- a/meshmode/discretization/poly_element.py +++ b/meshmode/discretization/poly_element.py @@ -40,6 +40,7 @@ Group types .. autoclass:: InterpolatoryQuadratureSimplexElementGroup .. autoclass:: QuadratureSimplexElementGroup .. autoclass:: PolynomialWarpAndBlendElementGroup +.. autoclass:: PolynomialRecursiveNodesElementGroup .. autoclass:: PolynomialEquidistantSimplexElementGroup .. autoclass:: LegendreGaussLobattoTensorProductElementGroup .. autoclass:: EquidistantTensorProductElementGroup @@ -50,6 +51,7 @@ Group factories .. autoclass:: InterpolatoryQuadratureSimplexGroupFactory .. autoclass:: QuadratureSimplexGroupFactory .. autoclass:: PolynomialWarpAndBlendGroupFactory +.. autoclass:: PolynomialRecursiveNodesGroupFactory .. autoclass:: PolynomialEquidistantSimplexGroupFactory .. autoclass:: LegendreGaussLobattoTensorProductGroupFactory """ @@ -208,12 +210,53 @@ class PolynomialWarpAndBlendElementGroup(_MassMatrixQuadratureElementGroup): polynomials in :math:`P^k`, hence usable for differentiation and interpolation. Interpolation nodes edge-clustered for avoidance of Runge phenomena. Nodes are present on the boundary of the simplex. + + Uses :func:`modepy.warp_and_blend_nodes`. + """ + @property + @memoize_method + def unit_nodes(self): + dim = self.mesh_el_group.dim + if self.order == 0: + result = mp.warp_and_blend_nodes(dim, 1) + result = np.mean(result, axis=1).reshape(-1, 1) + else: + result = mp.warp_and_blend_nodes(dim, self.order) + + dim2, nunit_nodes = result.shape + assert dim2 == dim + return result + + +class PolynomialRecursiveNodesElementGroup(_MassMatrixQuadratureElementGroup): + """Elemental discretization with a number of nodes matching the number of + polynomials in :math:`P^k`, hence usable for differentiation and + interpolation. Interpolation nodes edge-clustered for avoidance of Runge + phenomena. Nodes are present on the boundary of the simplex. + + Supports a choice of the base *family* of 1D nodes, see the documentation + of the *family* argument to :func:`recursivenodes.recursive_nodes`. + + Requires :mod:`recursivenodes` to be installed. + + .. [Isaac20] Tobin Isaac. Recursive, parameter-free, explicitly defined + interpolation nodes for simplices. + `Arxiv preprint `__. + + .. versionadded:: 2020.2 """ + def __init__(self, mesh_el_group, order, family, index): + super().__init__(mesh_el_group, order, index) + self.family = family + @property @memoize_method def unit_nodes(self): dim = self.mesh_el_group.dim - result = mp.warp_and_blend_nodes(dim, self.order) + + from recursivenodes import recursive_nodes + result = recursive_nodes(dim, self.order, self.family, + domain="biunit").T.copy() dim2, nunit_nodes = result.shape assert dim2 == dim @@ -342,6 +385,20 @@ class PolynomialWarpAndBlendGroupFactory(HomogeneousOrderBasedGroupFactory): group_class = PolynomialWarpAndBlendElementGroup +class PolynomialRecursiveNodesGroupFactory(HomogeneousOrderBasedGroupFactory): + def __init__(self, order, family): + self.order = order + self.family = family + + def __call__(self, mesh_el_group, index): + if not isinstance(mesh_el_group, _MeshSimplexElementGroup): + raise TypeError("only mesh element groups of type '%s' " + "are supported" % _MeshSimplexElementGroup.__name__) + + return PolynomialRecursiveNodesElementGroup( + mesh_el_group, self.order, self.family, index) + + class PolynomialEquidistantSimplexGroupFactory(HomogeneousOrderBasedGroupFactory): """ .. versionadded:: 2016.1 diff --git a/meshmode/mesh/generation.py b/meshmode/mesh/generation.py index c4bf260860d533ef39d7f188abd756aa215d97ac..c524292c1fffd741951b5838729014ac38db64b5 100644 --- a/meshmode/mesh/generation.py +++ b/meshmode/mesh/generation.py @@ -321,6 +321,9 @@ def make_group_from_vertices(vertices, vertex_indices, order, group_factory = SimplexElementGroup if issubclass(group_factory, SimplexElementGroup): + if order < 1: + raise ValueError("can't represent simplices with mesh order < 1") + el_origins = el_vertices[:, :, 0][:, :, np.newaxis] # ambient_dim, nelements, nspan_vectors spanning_vectors = ( diff --git a/requirements.txt b/requirements.txt index 818dd86c6ca8ad96a2ce29de55edbcc1f2fcdcba..574e50949645b761b6f3579c87fbb4e1dbfecc5e 100644 --- a/requirements.txt +++ b/requirements.txt @@ -1,21 +1,22 @@ numpy -git+https://gitlab.tiker.net/inducer/pytools.git -git+https://gitlab.tiker.net/inducer/gmsh_interop.git -git+https://gitlab.tiker.net/inducer/modepy.git -git+https://gitlab.tiker.net/inducer/pyopencl.git -git+https://gitlab.tiker.net/inducer/islpy.git +recursivenodes +git+https://gitlab.tiker.net/inducer/pytools.git#egg=pytools +git+https://gitlab.tiker.net/inducer/gmsh_interop.git#egg=gmsh_interop +git+https://gitlab.tiker.net/inducer/modepy.git#egg=modepy +git+https://gitlab.tiker.net/inducer/pyopencl.git#egg=pyopencl +git+https://gitlab.tiker.net/inducer/islpy.git#egg=islpy # required by pytential, which is in turn needed for some tests -git+https://gitlab.tiker.net/inducer/pymbolic.git +git+https://gitlab.tiker.net/inducer/pymbolic.git#egg=pymbolic # also depends on pymbolic, so should come after it -git+https://gitlab.tiker.net/inducer/loopy.git +git+https://gitlab.tiker.net/inducer/loopy.git#egg=loo.py # more pytential dependencies -git+https://gitlab.tiker.net/inducer/boxtree.git -git+https://gitlab.tiker.net/inducer/sumpy.git -git+https://gitlab.tiker.net/inducer/pytential.git@array-context +git+https://github.com/inducer/boxtree.git#egg=boxtree +git+https://github.com/inducer/sumpy.git#egg=sumpy +git+https://github.com/inducer/pytential.git#pytential # requires pymetis for tests for partition_mesh -git+https://gitlab.tiker.net/inducer/pymetis.git +git+https://gitlab.tiker.net/inducer/pymetis.git#egg=pymetis diff --git a/setup.py b/setup.py index 810d2560b674251f4b6b72d617e4fba74f38c692..1de30e57b763e1c01e91b5660e71023f93845238 100644 --- a/setup.py +++ b/setup.py @@ -44,6 +44,7 @@ def main(): "numpy", "modepy", "gmsh_interop", + "recursivenodes", "six", "pytools>=2020.3", "pytest>=2.3", diff --git a/test/test_meshmode.py b/test/test_meshmode.py index 9b1ec98739d55cb81529777bfaa79cc88e4d2440..5a77567f7dbe01bc468c887c7f700eb84568d015 100644 --- a/test/test_meshmode.py +++ b/test/test_meshmode.py @@ -22,6 +22,7 @@ OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. """ +from functools import partial from six.moves import range import numpy as np import numpy.linalg as la @@ -35,11 +36,12 @@ from meshmode.mesh import SimplexElementGroup, TensorProductElementGroup from meshmode.discretization.poly_element import ( InterpolatoryQuadratureSimplexGroupFactory, PolynomialWarpAndBlendGroupFactory, + PolynomialRecursiveNodesGroupFactory, PolynomialEquidistantSimplexGroupFactory, ) from meshmode.mesh import Mesh, BTAG_ALL from meshmode.array_context import PyOpenCLArrayContext -from meshmode.dof_array import thaw, flat_norm, flatten +from meshmode.dof_array import thaw, flat_norm, flatten, unflatten from meshmode.discretization.connection import \ FACE_RESTR_ALL, FACE_RESTR_INTERIOR import meshmode.mesh.generation as mgen @@ -207,7 +209,9 @@ def test_box_boundary_tags(dim, nelem, group_factory): @pytest.mark.parametrize("group_factory", [ InterpolatoryQuadratureSimplexGroupFactory, - PolynomialWarpAndBlendGroupFactory + PolynomialWarpAndBlendGroupFactory, + partial(PolynomialRecursiveNodesGroupFactory, family="lgl"), + #partial(PolynomialRecursiveNodesGroupFactory, family="gc"), ]) @pytest.mark.parametrize("boundary_tag", [ BTAG_ALL, @@ -300,7 +304,7 @@ def test_boundary_interpolation(ctx_factory, group_factory, boundary_tag, print(eoc_rec) assert ( eoc_rec.order_estimate() >= order-order_slack - or eoc_rec.max_error() < 1e-14) + or eoc_rec.max_error() < 3e-14) # }}} @@ -1442,6 +1446,39 @@ def test_mesh_multiple_groups(ctx_factory, ambient_dim, visualize=False): assert error < 1.0e-11, error +def test_array_context_np_workalike(ctx_factory): + ctx = ctx_factory() + queue = cl.CommandQueue(ctx) + actx = PyOpenCLArrayContext(queue) + + from meshmode.mesh.generation import generate_regular_rect_mesh + mesh = generate_regular_rect_mesh( + a=(-0.5,)*2, b=(0.5,)*2, n=(8,)*2, order=3) + + from meshmode.discretization import Discretization + from meshmode.discretization.poly_element import \ + PolynomialWarpAndBlendGroupFactory as GroupFactory + discr = Discretization(actx, mesh, GroupFactory(3)) + + for sym_name, n_args in [ + ("sin", 1), + ("exp", 1), + ("arctan2", 2), + ("minimum", 2), + ("maximum", 2), + ("where", 3), + ]: + args = [np.random.randn(discr.ndofs) for i in range(n_args)] + ref_result = getattr(np, sym_name)(*args) + + actx_args = [unflatten(actx, discr, actx.from_numpy(arg)) for arg in args] + + actx_result = actx.to_numpy( + flatten(getattr(actx.np, sym_name)(*actx_args))) + + assert np.allclose(actx_result, ref_result) + + if __name__ == "__main__": import sys if len(sys.argv) > 1: